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Abstract 

Although  the  theory  of  linear  control  systems  is  highly  mature,  nonlinear  control-system  design 
techniques  remain  relatively  undeveloped.  In  real-world  applications  such  as  vibration  suppression 
for  flexible  structures  and  large  angle  rigid-body  spacecraft  maneuvers,  nonlinear  plants  generally 
require  nonlinear  controllers,  while  linear  plants  often  benefit  from  the  implementation  of  nonlinear 
controllers  in  the  presence  of  structured  plant  uncertainty,  actuator  constraints,  and  nonquadratic 
performance  criteria.  This  report  discusses  progress  in  several  areas  relating  to  the  role  of  non- 
linearities  in  feedback  control.  These  areas  include  Lyapunov  function  theory,  chaotic  controllers, 
Statistical  Energy  Analysis,  phase  robustness,  and  optimal  nonlinear  feedback  control. 
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1.0  Introduction 


1.1  Review  of  Linear  Multivariable  Control  Theory 

The  bulk  of  current  control-system  practice  is  based  upon  linear  control  theory.  Classical  single¬ 
loop  design  methods,  whose  basic  development  predates  1960,  are  widely  utilized  in  practice.  For 
high-performance  multi-loop  applications,  modem  multivariable  techniques  are  finding  their  way 
into  practice.  A  broad  spectrum  of  linear  multivariable  control  techniques  has  reached  the  graduate 
curriculum  (see,  for  example,  the  in-depth  textbook  [1]).  Moreover,  the  advanced  development  of 
such  methods  is  reflected  in  the  availability  of  several  computer-aided  design  packages. 

Within  modem  multivariable  control  theory  there  are  severed  major  thrusts  of  development 
that  can  be  identified.  From  a  state  space  perspective,  the  original  work  of  Kalman  and  others 
has  led  to  a  rather  complete  theory  of  Hj-optimal  linear-quadratic- Gaussian  (LQG)  control  design 
[2-4].  Furthermore,  an  elegant  state  space  theory  within  a  geometric  rather  than  optimization 
framework  has  been  developed  in  [5].  Multivariable  extensions  of  classical  frequency-domain  ideas 
have  undergone  significant  development  along  a  number  of  paths.  For  example,  classical  ideas  have 
been  generalized  to  the  multivariable  setting  in  [6,7],  while  an  optimal  design  theory  based  upon  a 
frequency-domain  (Hoo)  criterion  was  pioneered  in  [8]  and  further  developed  in  numerous  papers 
(see,  e.g.,  [9,10]and  references  therein).  We  also  note  the  development  of  further  sophisticated 
approaches  within  an  algebraic  transfer  function  setting  [11,12]. 

It  is  also  worthwhile  reviewing  some  recent  trends  in  linear  multivariable  control,  namely,  robust 
control  and  controller  simplification.  Robust  control  refers  to  the  need  to  effect  desired  closed-loop 
performance  (e.g.,  tracking  and  disturbance  rejection)  in  spite  of  plant  modeling  uncertainties. 
Within  classical  theory,  the  related  concept  of  sensitivity  plays  a  key  role,  while  multivariable 
problems  require  more  sophisticated  approaches.  Numerous  robust  control-design  techniques  have 
been  developed  under  a  variety  of  assumptions  concerning  the  plant  uncertainty.  Unstructured 
uncertainty  is  addressable  via  Hoo  methods  [9,10],  while  specialized  techniques  are  required  for 
more  highly  structured  plant  uncertainty;  see,  e.g.,  [13-16,1.20,1.30].  In  addition,  recent  results 
concerning  the  state  space  solution  of  Hoo  problems  yield  greater  unification  of  state  space  and 
frequency  domain  synthesis  techniques  [17-19,1.29] . 

The  second  trend  in  linear  multivariable  control  theory  we  note  here  involves  controller  sim¬ 
plification  issues.  While  modem  design  techniques  such  as  LQG  theory  produce  high-order  con- 
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trailers,  it  is  desirable  in  practice  to  employ  the  simplest  controller  meeting  design  specifications. 
Here,  “simplicity”  may  refer  to  dynamic  dimension,  number  of  digital  operations,  degree  of  de¬ 
centralization,  and  other  considerations  affecting  implementation,  cost,  reliability,  etc.  References 
[20-24,1.23,11.89,11.91]  are  representative  of  progress  made  in  this  area. 

1.2  Nonlinear  Control  Theory 

“Nonlinear  control  theory”  refers  to  control  theory  in  which  either  the  controller  or  the  plant  (or 
both)  is  nonlinear.  This  theory  is  not  as  extensively  developed  as  linear  multivariable  control  the¬ 
ory.  The  principal  approaches  to  nonlinear  multivariable  control  design  include  local  linearization, 
global  linearization,  the  second  method  of  Lyapunov,  variable  structure  control,  optimization-baaed 
methods,  and  differential-geometric  methods.  With  these  general  classifications  in  mind,  we  can 
identify  several  advantages  of  nonlinear  control  over  linear  control.  In  this  regard  it  is  useful  to 
consider  three  cases  in  which  the  theory  is  applied  (see  Figure  1.2-1):  (i)  nonlinear  control  for  linear 
plants,  (n)  linear  control  for  nonlinear  plants,  and  (tit)  nonlinear  control  for  nonlinear  plants. 

The  role  of  nonlinearities  in  control  theory  can  best  be  understood  by  reviewing  the  assumptions 
and  limitations  of  standard  linear-quadratic-Gaussian  (LQG)  theory.  As  its  name  implies,  LQG 
theory  is  based  upon  three  fundamental  assumptions  (Figure  1.2-2) 

•  the  plant  dynamics  and  measurement  equations  are  linear  in  both  the  state  and  control 
variables 

•  the  performance  measure  to  be  minimized  is  quadratic 

•  the  plant  disturbances  and  measurement  noise  are  additive  Gaussian  white  noise 
In  addition  to  these  explicit  assumptions  the  following  implicit  assumptions  are  crucial: 

•  the  plant  model  is  completely  accurate 

•  mean-square  control  effort  is  limited 

Under  these  assumptions,  a  major  result  of  modern  control  theory  [2]  states  that  the  opti¬ 
mal  controller  is  given  by  the  linear  controller  consisting  of  the  Wiener-Kalman  filter  followed  by 
the  optimal  linear-quadratic  regulator.  Hence  in  this  case  nonlinear  controllers  cannot  improve 
performance. 
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•  EXTENDED  DOMAIN  OF  ATTRACTION 

•  NONQUADRATIC  COST 


Figure  1.2—1.  In  nonlinear  control  theory,  the  plant  and/or  controller  is  nonlinear. 
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•  QUADRATIC  COST 

•  LINEAR  DYNAMICS  AND  MEASUREMENTS 

•  ADDITIVE  GAUSSIAN  DISTURBANCES 

•  NO  MODELING  UNCERTAINTY 

•  MEAN-SQUARE  ACTUATOR  BOUNDS 

* 

•  LINEAR  CONTROLLER  IS  OPTIMAL 
(LOG  THEORY) 


•  NONLINEAR  DYNAMICS  AND/OR  MEASUREMENTS 

•  NONGAUSSIAN,  NONADDITIVE  DISTURBANCES 

•  MODELING  UNCERTAINTY 

•  NONQUADRATIC  COST 

•  AMPLITUDE  ACTUATOR  BOUNDS 


•  NONLINEAR  CONTROLLER  IS  OPTIMAL 


Figure  1.2—2.  Linear  controllers  are  generally  optimal  for  only  a  narrow  class  of  linear-quadratic- 
Gaussian  problems. 
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Suppose,  however,  that  not  all  of  the  assumptions  of  LQG  theory  are  valid  for  a  given  problem, 
that  is,  one  or  more  of  the  following  conditions  applies: 

•  the  plant  dynamics  and/or  measurement  equation  is  nonlinear 

•  the  disturbances  are  either  nonadditive  or  non-Gaussian 

•  the  relevant  performance  measure  is  nonquadratic 

•  the  plant  model  is  uncertain 

•  control  effort  is  limited  by  amplitude  (Loo)  or  total  fuel  (Li)  constraints. 

In  real  world  applications,  of  course,  all  of  these  conditions  apply,  at  least  to  some  extent.  The 
actual  extent  to  which  each  one  must  be  considered  is  problem-dependent.  In  each  of  these  cases 
there  is  no  reason  to  expect  that  a  linear  controller  is  optimal  or  even  appropriate.  Nevertheless,  it  is 
still  desirable  for  a  variety  of  reasons  to  seek  linear  controllers,  and  much  of  control  theory  has  been 
directed  toward  this  goal.  Ultimately,  however,  we  are  faced  with  the  following  question:  When 
is  it  necessary  or  advantageous  to  implement  nonlinear  controllers  in  place  of  linear  controllers? 
Nonlinear  controllers  will  generally  entail  more  difficult  performance  validation  and  implementation 
complexity  (Figure  1.2-3).  Furthermore,  we  note  that  an  additional  level  of  controller  complexity 
involves  time- varying  control  (for  either  linear  or  nonlinear  controllers)  (Figure  1.2-4).  We  now 
examine  the  possible  benefits  of  nonlinear  controllers. 

1.3  Linear  Versus  Nonlinear  Controllers 

Let  us  first  consider  the  problem  of  nonlinear  plant  dynamics.  Such  nonlinearities  arise  in  a 
wide  variety  of  engineering  applications  [25],  Nevertheless,  linear  control  theory  has  been  developed 
to  deal  with  large  classes  of  nonlinearities,  for  example,  as  bounded  by  a  sector  [26-28].  Lure’s 
problem,  the  Aizermann  conjecture,  and  the  circle  and  Popov  criteria  are  all  traditional  control 
theory  topics  dealing  with  nonlinearities. 

In  many  applications,  however,  the  nonlinearities  are  well  modeled  to  the  extent  that  their 
detailed  structure  can  be  exploited  in  control  design.  For  example,  in  the  case  of  a  single  rigid 
body  we  have  Euler’s  equation 

Jd)  +  <jJ  x  Ju>  =  /(u), 

where  J  denotes  the  moment  of  inertia,  w  denotes  angular  velocity,  and  /( u)  denotes  applied  torque. 
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*  GREATER  GENERALITY 

*  IMPROVED  PERFORMANCE 
FOR  MODELING  ACCURACY 

*  MORE  COMPLEX  TO 
IMPLEMENT 

*  HARDER  TO  VALIDATE 


► 


LINEAR 

TIME-INVARIANT 

CONTROLLER 

GAIN  SCHEDULED 
LTI  CONTROLLER 

NONLINEAR 

CONTROLLER 

ADAPTIVE 

CONTROLLER 

◄ 


*  MORE  RESTRICTIVE  CONTROLLER  CLASS 

*  PERFORMANCE  LIMITED  BY  MODELING  ACCURACY 

*  SIMPLER  TO  IMPLEMENT 

*  EASIER  TO  VALIDATE 


DESIGN  GUIDELINES 

*  TRY  TO  MEET  PERFORMANCE  SPECIFICATIONS 
WITH  SIMPLEST  POSSIBLE  CONTROL  LAW 

*  IF  SPECIFICATIONS  CANNOT  BE  MET,  THEN 
INCREASE  CONTROL  LAW  COMPLEXITY  AND 
ASSESS  PERFORMANCE/IMPLEMENTATION/ 
VALIDATION  TRADEOFFS 


Figure  1.2-3.  Nonlinear  controllers  offer  improved  performance,  but  may  entail 
greater  implementation  and  validation  complexity. 
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Figure  1.2-4.  Linear  and  nonlinear  controllers  may  be  either  time-invariant  or  time- 
varying. 
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The  quadratic  gyroscopic  term  ui  x  Joj  is  significant  in  rapid  maneuvers  involving  large  structures. 
Since  the  structure  of  the  nonlinearity  in  this  case  is  crucial,  we  expect  nonlinear  controllers  to  play 
a  role  [29-57]. 

Next  consider  the  problem  of  optimizing  a  nonquadratic  performance  measure.  In  this  case  it 
can  generally  be  expected  that  linear  controllers  are  not  optimal.  Time-optimal  performance  leads 
to  bang-bang  controllers,  which  are  nonlinear,  while  higher-order  (polynomial)  performance  mea¬ 
sures  lead  to  higher-order  feedback  laws  [58-68] .  For  example,  consider  the  effect  of  a  nonquadratic 
performance  measure  as  addressed  in  [68].  As  shown  in  Figure  1.2-5,  a  super-linear  state  feedback 
(in  this  case  a  quadratic  control)  can  efficiently  regulate  small  amplitude  signals,  even  driving  the 
state  to  zero  in  finite  time  (if  one  neglects  measurement  and  disturbance  noise  effects).  A  theory 
of  sublinear  control  for  finite-time  control  is  developed  in  [69].  An  additional  performance  aspect 
is  transient  behavior  [70]  which  is  difficult  to  capture  by  means  of  scalar  performance  measures. 

There  are,  however,  nonquadratic  performance  measures  for  which  the  optimal  controller  ia 
linear.  In  particular,  this  is  the  case  for  Hoo  optimal  control.  For  this  problem  the  goal  is  to 
minimize  the  worst-case  disturbance  attenuation  over  all  frequencies.  The  Hoo  problem  differs 
mathematically  from  the  LQG  problem  due  to  the  modeling  of  disturbances  and  error  signals  as 
deterministic  La  functions.  Connections  with  the  LQG  setting  can  be  established  by  means  of  an 
exponential-of-quadratic  performance  functional  with  white  noise  disturbances  [71-84], 

Problems  involving  uncertain  plant  models  have  motivated  the  subject  of  robust  control  theory. 
One  approach  to  robust  control  involves  modeling  the  uncertainty  by  means  of  the  Hoo  norm  and 
then  applying  theory  to  guarantee  robust  stability  and  performance.  In  this  case  and  for 
related  problems  in  robust  control,  it  has  been  shown  that  nonlinear  controllers  offer  no  advantage 
over  linear  controllers  [85-90].  Though  valuable,  these  results  consider  only  restricted  uncertainty 
characterizations  (e.g.,  unstructured  uncertainty)  [86-88],  very  special  performance  measures  (e.g., 
Hoo  performance)  [85],  or  limited  definitions  of  stability  (e.g.,  quadratic  stability)  [90].  In  fact,  from 
the  previous  discussion  on  the  optimality  of  nonlinear  controllers  for  nonquadratic  performance 
criteria,  it  is  reasonable  to  conjecture  that  for  a  variety  of  system  performance  measures  nonlinear 
controllers  can  yield  better  robust  performance  than  linear  controllers.  In  fact,  it  is  even  possible 
that  the  controller  that  solves  the  robust  quadratic  performance  problem 

roo 

max  /  (xtRiX  +  uTR2u)dt 
(AA,AB)eUj0 

x(t)  =  (A  +  AA)x{t)  +  (B  +  dfl)tt(t), 


min< 

«(‘)  l 
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p=|  :  u  =  -2|]x||x  =*■  EFFICIENT  REGULATION  FOR  LARGE  AMPLITUDES 

p  =  1  :  u  =  -2x  ==>  EXPONENTIAL  DECAY  (LINEAR) 

:  u=  -2||x|r^x=>  EFFICIENT  REGULATION  FOR  SMALL  AMPLITUDES 

4 

Figure  1.2-5.  As  shown  in  [68],  nonlinear  controls  can  be  shaped  to  give  efficient  regulation 
for  various,  selected  vibration  amplitude  regimes. 
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is  a  nonlinear  controller.  Results  that  indicate  that  nonlinear  controllers  can  yield  improved  ro¬ 
bustness  properties  are  given  in  [91-94]. 

An  additional  advantage  of  nonlinear  controllers  is  the  ability  to  address  actuator  saturation 
limitations.  In  practice  any  electromechanical  device  used  as  a  control  actuator  is  subject  to 
limitations  on  maximum  force,  torque  output,  power  consumption,  stroke,  and  angular  speed  limits. 
Thus,  in  reality,  control-design  optimization  must  account  for  constraints  on  the  maximum  value 
of  actuator  force  or  similar  constraints  on  internal  signals  associated  with  the  actuator  dynamics. 
The  simplest  such  realistic  constraint  takes  the  form  of  a  pointwise  bound  on  the  actuator  force 
output,  i.e., 

MOI  ^  «msx, 

where  umkX  is  the  largest  physically  possible  magnitude  of  the  actuator  output.  The  above  pointwise 
bound  is  an  example  of  an  design  constraint  and  differs  crucially  from  L2  constraints  such  as 

IE[us]  <  a * 


or 

f  u <  a a, 

J  o 

which  correspond  to  power  and  energy  constraints,  respectively. 

Figure  1.2-6  illustrates  the  ramifications  of  pointwise  bounded  actuator  constraints.  Suppose 
that  the  plant  is  linear  except  for  the  physical  constraint  |u(i)|  <  umax  and  that  the  system  is 
subject  to  an  initial  impulse  disturbance.  If  one  designs  an  optimal  regulator  using  the  integral 
square  condition  (for  analytical  convenience)  as  the  constraint  on  the  optimization  problem,  then 
the  resulting  controller  is  linear.  Moreover,  one  can  choose  linear  gains  such  that  the  peak  actuator 
output  is  less  than  the  physically  imposed  limit  of  umax.  Then  (see  the  top  half  of  Figure  1.2-6), 
following  the  initial  disturbance,  all  signals,  including  the  actuator  output,  decay  exponentially. 
Note  that  although  u(t)  just  satisfies  the  physical  constraint  [u(t)|  <  umax  for  small  t,  for  larger 
t,  |u(t)|  is  small.  Thus,  in  this  case  actuator  capability  is  wasted.  On  the  other  hand  (see  bottom 
half  of  Figure  1.2-6),  it  is  possible  to  design  a  sublinear  feedback  control  for  which  the  actuator 
uses  nearly  its  full  capacity  while  the  system  state  is  driven  to  zero  faster  than  exponentially.  In 
fact,  it  is  well  known  that  minimal-time  maneuvers  actually  require  bang-bang  control.  Variable 
structure  (nonlinear)  controllers,  which  can  be  viewed  as  generalizations  of  bang-bang  controllers, 
can  also  be  used  to  control  linear  systems  while  efficiently  utilizing  actuator  capabilities. 
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UNEAR  OPTIMAL  SATISFYING  J^Vdt  Su2^ 


FOR  SMALL  SIGNALS,  U(t) 
SMALL  — ►  ACTUATOR 
CAPABILITY  WASTED! 

"SUBOPT1MAL”  NONUNEAR  CONTROL 


t* 


FINITE  TIME 
TO  ZERO 


*  NONUNEAR  CONTROLS  CAN  MORE  EFFICIENTLY 
*  *  UTILIZE  REALISTIC  ACTUATOR  CAPABILITIES 


Figure  1.2-6.  Nonlinear  controllers  can  utilize  actuators  more  efficiently  than  linear  con* 
trollers  in  the  presence  of  saturation  bounds. 
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A  specialized  class  of  nonlinear  controllers  for  linear  plants  is  the  class  of  adaptive  controllers. 
In  contrast  to  fixed-gain  controllers,  which  maintain  prespecified  constants  within  the  feedback  law, 
adaptive  controllers  adjust  feedback  gains  to  improve  closed-loop  stability  and  performance  when 
the  plant  is  uncertain.  Adaptive  controllers  generally  utilize  probing  signals  to  excite  the  plant 
dynamics  and  thereby  identify  plant  parameters.  Feedback  gains  can  then  be  adjusted  to  account 
for  the  identification  data.  The  overall  process  of  identification  and  adjustment  clearly  constitutes  a 
nonlinear  control  law.  Thus,  the  adaptive  control  literature  can  be  viewed  as  a  specialized  subclass 
of  nonlinear  control,  although  for  historical  reasons  this  categorization  is  rarely  utilized.  For  our 
purposes,  viewing  adaptive  controllers  as  nonlinear  controllers  is  particularly  useful.  For  example, 
as  discussed  above,  nonlinear  controllers  can  be  viewed  as  a  specialized  form  of  robust  controllers 
for  uncertain  linear  plants. 

The  distinction  between  nonlinear  controllers  and  adaptive  controllers  has  narrowed  in  recent 
years  with  the  development  of  adaptive  controllers  not  requiring  explicit  probing  signals  [95-102]. 
These  results  show  that  there  exist  nonlinear  controllers  that  can  stabilize  generic  classes  of  sys¬ 
tems  characterized  by  minimal  a  priori  data.  Although  these  controllers  are  usually  thought  of  as 
adaptive  since  the  feedback  gains  are  continually  adjusted,  the  feedback  laws  are  clearly  nonlinear 
controllers  of  special  structure. 

1.4  Overview  of  this  Report 

The  central  result  of  control  system  analysis  and  design  is  Lyapunov’s  method.  The  ability 
to  construct  a  positive-definite  functional  that  decays  aloqg  system  trajectories  is  sufficient  to 
guarantee  asymptotic  stability.  Design  via  Lyapunov  functions  need  not  be  associated  with  the 
optimization  of  a  performance  measure  although,  as  discussed  in  Section  5,  the  converse  is  often 
true,  that  is,  optimal  design  may  be  predicated  on  a  Lyapunov  function.  Hence,  in  our  view, 
Lyapunov’s  method  ultimately  comprises  the  most  fundamental  technique  in  nonlinear  (as  well  as 
linear)  control  theory. 

This  program  is  thus  focussing  on  several  problem  areas  relating  to  Lyapunov  theory.  The 
interrelationships  among  these  areas  is  shown  in  Figure  1.4-1.  In  Section  2  we  describe  progress 
in  analyzing  energy  flow  in  coupled  mechanical  systems.  The  results  obtained  thus  far  extend  the 
foundations  of  Statistical  Energy  Analysis.  In  Section  3  we  apply  the  results  of  Section  2  along 
with  applications  to  the  design  of  chaotic  controllers  for  enhanced  energy  dissipation.  Section  4  is 
devoted  to  progress  in  developing  a  theory  of  robustness  due  to  phase  properties.  Finally,  Section 
5  discusses  optimal  nonlinear  control  theory. 
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Figure  1.4-1.  The  program  focusses  on  several  related  research  problems  relevant  to  nonlinear 
control. 
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2.0  Power  Flow  and  Statistical  Energy  Analysis 


2.0  Energy  Flow  and  Statistical  Energy  Analysis 


It  is  well  known  from  thermodynamics  that  energy  flows  from  hot  objects  to  cold  objects. 
It  is  less  well  known,  however,  that  a  similar  phenomenon  occurs  in  coupled  mechanical  systems 
with  modal  energy  playing  the  role  of  temperature.  Energy  and  power  flow  concepts,  often  called 
Statistical  Energy  Analysis  (SEA),  have  proven  to  be  useful  tools  for  analyzing  linear  dynamic 
systems  [237-249].  Hence  this  phase  of  the  program  is  devoted  to  the  further  development  of  these 
ideas  to  support  nonlinear  analysis  and  design.  In  Section  3  these  ideas  are  used  to  analyze  and 
design  chaotic  feedback  controllers. 

2.1  Energy  Flow  in  Coupled  Dynamical  Systems 

The  objective  of  SEA  is  to  model  energy  flow  among  coupled  dynamical  subsystems.  SEA 
was  originally  developed  for  acoustical  analysis  involving  very  large  numbers  of  modes  that  may 
be  poorly  modeled.  Many  of  the  concepts  of  SEA  as  applied  to  high  dimensional  systems  (such  as 
equipartition  of  energy)  have  close  connections  with  statistical  mechanics  of  many  particle  systems. 
Although  SEA  theory  has  been  widely  applied,  rigorous  analytical  results  have  been  available  only 
for  identical  couplings  or  for  weak  interactions.  Under  this  program  we  have  extended  SEA  theory 
to  address  an  arbitrary  number  of  subsystems  with  arbitrary  coupling. 

In  this  section  we  summarize  results  on  SEA  which  are  developed  in  the  paper  entitled  “Power 
Flow,  Energy  Balance,  and  Statistical  Energy  Analysis  for  Large-Scale  Interconnected  Systems.” 
This  paper,  which  contains  all  details  of  the  results  reported  here,  appears  in  Appendix  C. 

To  summarize  these  results  consider  the  system 

x  =  Ax  +  Gx  +•  w,  (1) 

where  the  state  z  €  Cn,  the  uncoupled  dynamics  matrix  A  is  given  by 

A  =  - u  +  jf]+H , 

v  =  diag  (yit  ...,*/„)€  IRnXn,  */<  >  0, 

H=  diag  (Hi,.  .,#„)€  <D"X", 
diag(tfi,...,flre)eIR"x", 
and  where  G  denotes  the  coupling  among  subsystems,  that  is, 

G  €  <DnXn,  Ga  =  0,  i=l,...,n. 
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The  additive  forcing  w(t)  is  taken  to  be  white  noise  with  intensity  V  >  0. 

The  first  step  of  our  approach  is  to  note  that  for  an  output  signal 

z  —  Cx,  (2) 

the  steady-state  mean-square  response  is  given  by 

J  =  lim  IE[z*z] 

t-°°  (3) 

=  tr[CTCQ], 

where 

Q=  lim  IEfix*]. 

t— »oo  1  1 

It  is  well  known  that  the  steady-state  covariance  Q  is  given  by  the  algebraic  Lyapunov  equation 


0  =  AQ  +  QA*+GQ  +  QG*  +  V. 


ha  many  practical  situations  it  can  be  argued  (see  Appendix  C)  that  the  principal  contribution 
to  J  is  due  to  the  diagonal  elements  of  Q.  Hence  our  main  result  is  based  on  a  direct  characterization 
of  the  diagonal  elements  of  Q  in  terms  of  V,  which  is  obtained  by  eliminating  the  off-diagonal 
elements  of  Q.  To  do  this,  we  rewrite  (4)  as 

0  =  A{Q}  +  { Q}A •  +  {G(Q»  +  {(Q)G*>  +  {V},  (5) 

0  =  A(Q)  +  {( Q)A *  +  {G{Q))  +  «Q)G*>  +  G{Q}  +  (Q}G*,  (6) 


where  {•}  and  (•)  denote  the  diagonal  part  and  off-diagonal  part  of  a  matrix,  respectively.  Here  we 
have  assumed  for  convenience  that  (V)  =  0. 

Next,  we  apply  Kronecker  matrix  algebra  to  solve  (5),  (6)  for  {Q}  in  terms  of  {V}.  To  state 
the  main  result  define  the  vector  E  of  steady-state  mean-square  state  energy's 


where  ^  =  Qu,  i  =  1, . . . , n,  and  the  vector 
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which  corresponds  to  {V}.  Then  we  obtain  the  following  consequence  of  (5),  (6): 


o 1  +  P)E  =  V, 


(V 


where 

H  =  diag  {2vi  -  2Re(l?i), . . . ,  2i/n  -  2R e(Hn)},  (8) 

and 

P  =  iT{&  ©  G)fx[(A  ©  A)  ©  (G  ©  G))~l£±{G  ©  G)£.  (9) 

In  (8)  and  (9),  “Re”  denotes  real  part,  ©  denotes  Kronecker  sum,  and  £,  £±,  and  £  denote  n3  x  n3 
matrices  of  special  structure  whose  element  are  ones  and  zeros.  It  can  be  shown  that  P  is  real. 


Tb  elucidate  the  meaning  of  (7)  we  can  write  its  fcth  component  as 


(2 Vk  -  2Re  Hk)Ek 

power  dissipated 
by  the  Arth  mode 
due  to  damping 


+  nk 

power  flow 
from  the  fcth  mode 
to  all  other  modes 
due  to  coupling 


Vfcfc, 


power  input 
due  to  external 
disturbances 


(10) 


where  Ilk  has  the  form 

n 

nh  =  YtpkiEi,  Pki  €  nt.  (ii) 

t=i 

The  matrix  P  can  be  viewed  as  a  power  flow  matrix,  while  relation  (10)  thus  has  the  form  of  a 
power  flow  equation.  To  arrive  at  an  energy  balance  relation  we  consider  the  case  in  which 


Pu<0,  k^t,  k,t=  l,...,n. 


(12) 


This  occurs,  for  example,  if  the  subsystem  coupling  is  sufficiently  weak.  If,  in  addition,  the  couplings 
are  energy  conservative  (for  example,  passive),  then  it  can  be  shown  that 


Pkk  =  '£,\Pki\,  k=  1,. . .  ,n.  (13) 

tml 

l+h 


Then,  defining  out  =  |  Piu\,  k  ^  l,  so  that  <ru  >  0,  it  follows  from  (11)  that 


Hie  =  ^,Okt{Ek  ~  Ei).  (14) 

<-L 

t+h 

In  other  words,  power  flow  from  the  fcth  mode  to  all  other  modes  is  the  sum  of  the  individual 
power  flows  from  mode  k  to  mode  t,  which  are  proportional  to  the  energy  differences  Ek  -  Et . 
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Note  that  power  always  flows  from  more  energetic  modes  to  less  energetic  modes  (because  of  the 
nonnegativity  of  the  coefficients  Ckt)-  Substituting  (14)  into  (7)  yields 

n 

HkEk  +  ^  Okt{Ek  -  Et)  =  Vfe,  (15) 

t+k 

which  is  an  energy  balance  relation.  Equations  (10)  and  (15),  which  govern  energy  exchange  among 
coupled  oscillators,  are  completely  analogous  to  the  equations  of  thermal  transfer  with  the  modal 
energies  playing  the  role  of  temperatures. 

In  physical  situations  involving  nonconservative  couplings,  we  have  shown  that  although  (13) 
no  longer  holds,  it  is  still  possible  in  the  case  of  weak  couplings  to  obtain  a  generalized  power  flow 
proportionality.  In  this  case  there  exists  a  set  of  positive  scale  factors  Dk  >  0,  k  =  1, . . . ,  n,  such 
that,  with  Eh  =  jj^Ei,,  the  energy  difference  power  flow  proportionality  is  given  by 

n 

=  5>w(Efc  -  &),  (16) 

£—  1 

t+h 

where  &kl  —  Di<ru .  Note  that  (16)  is  not  merely  a  rewriting  of  (14)  since  in  general  Dk  /  Dt. 
With  (16),  the  energy  equation  (7)  assumes  the  form  of  a  generalized  energy  balance  relation  given 

by 

fo&k  +  £  *kl{Ek  -  Et)  =  Vk,  (17) 

im  I 
t+k 

where  k  —  1, . . . ,  n.  That  is,  there  is  a  set  of  re-scaled  energies  such  that  (7)  looks  like  the  equations 
of  thermal  transfer. 

Furthermore,  while  deriving  energy  difference  power  flow  proportionality  relations,  we  have 
also  shown  that  the  explicit  expressions  given  for  the  power  flow  matrix  P  in  the  SEA  literature 
are  actually  first-term  approximations  in  a  series  expansion  for  P.  Indeed,  it  turns  out  that  P, 
which  is  given  by  a  complicated  expression  involving  v,n,H,  and  G,  agrees  with  the  customary 
SEA  expressions  for  "small”  G.  This  in  done  by  obtaining  explicit  expressions  for  the  terms  of  a 
series  expansion  of  P  in  ascending  powers  of  the  matrix  elements  of  G. 

Since  the  modal  energies  obey  equations  analogous  to  those  of  thermal  transfer,  it  might  be 
expected  that  if  the  coupling  coefficients  Gki  are  large  compared  to  the  modal  dampings,  then  the 
energies  should  be  approximately  equal,  that  is, 

Ei  —  Ei  ~  ~  En.  (18) 
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The  paper  in  Appendix  C  provides  a  formulation  and  proof  of  this  “energy  equipartitioning”  phe¬ 
nomenon. 
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3.0  Chaotic  Controllers 


3.0  Chaotic  Controllers 


In  Section  2  we  explored  the  notions  of  power  flow  and  energy  balance  in  interconnected  systems. 
Our  next  goal  is  to  apply  these  ideas  to  the  analysis  and  design  of  nonlinear  feedback  controllers. 

lb  do  this  we  need  only  view  the  plant  and  controller  as  a  pair  of  interacting  subsystems.  If 
disturbance  rejection  is  an  objective,  then  we  seek  to  design  a  controller  that  maximizes  power  flow 
from  the  plant  to  the  controller.  Within  an  Hoo  context  this  idea  has  been  explored  in  the  recent 

paper 

D.  MacMartin  and  S.  R.  Hall,  “An  Hoo  Power  Flow  Approach  to  Control  of  Uncertain  Struc¬ 
tures,”  Proe.  Amer.  Contr.  Conf.,  pp.  3073-3080,  San  Diego,  CA,  May  1990. 

One  of  the  main  ideas  discussed  in  this  paper  is  that  power  flow  out  of  the  structure  is  maximized 
to  the  extent  that  the  controller  is  able  to  match  the  impedance  of  the  plant. 

In  this  section  we  develop  a  nonlinear  controller  that  exploits  the  phenomenon  of  chaos.  Our 
principal  goal  is  to  demonstrate  that  this  controller  can  enhance  power  flow  from  the  plant  to  the 
controller  by  introducing  nonlinearities  that  induce  broadband  spectral  properties  in  the  controller. 
A  power  flow  analysis  is  then  used  to  show  that  energy  can  be  transferred  more  efficiently  between 
arbitrary  plant  and  compensator  modes.  Details  of  these  results  are  given  in  “A  Nonlinear  Vibration 
Control  Design  with  a  Neural  Network  Realization”  which  appears  in  Appendix  D. 

3.1  Turbulence  Model  for  Chaotic  Controller  Design 

A  unique  feature  of  nonlinear  systems  is  the  energy  cascade  mechanism  illustrated  in  Figure 
3.1-1.  Here,  energy  originally  injected  within  some  lower  frequency  band  can  be  dispersed  to  higher 
frequency  bands  by  virtue  of  coupling  among  the  vibrational  modes  of  the  structure.  Eventually 
the  energy  is  transferred  to  very  high  frequencies  where  it  is  dissipated  into  heat  by  means  of 
natural  structural  damping  or  by  the  action  of  an  additional  energy  dissipative  control  law.  Thus, 
a  nonlinear  controller  such  as  illustrated  in  Figure  3.1-1  can  be  viewed  as  a  catalyst  for  transmuting 
vibration  more  rapidly  into  heat. 

The  controller  illustrated  in  Figure  3.1-1  can  be  realized  by  a  purely  mechanical  device  con¬ 
sisting  of  a  chamber  containing  a  number  of  particles  of  given  mass  that  undergo  free  translational 
motion  except  for  collisions  with  the  chamber  walls  and  with  one  another.  This  is  essentially  the 
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ENERGY  CASCADES  ENERGY  EVENTUALLY 

INTO  HIGHER  DISSIPATES  INTO 

FREQUENCY  BANOS  HEAT  AT  HIGH  FREQUENCY 


•  ENERGY  CASCADE  MECHANISM  OFFERS  POTENTIAL  FOR  EXTREMELY 
ROBUST,  RAPID  ATTENUATION  OF  LOWER  FREQUENCY  VIBRATION. 


•  SUCH  NONLINEAR  CONTROLLERS  ARE  A  CATALYST  FOR  TRANSMUTING 
VIBRATION  MORE  RAPIDLY  INTO  HEAT. 


•  COMPENSATOR  (BY  ITSELF)  COULD  BE  CHAOTIC.  BUT  WHEN 
INTERCONNECTED  WITH  THE  PLANT,  ITS  DAMPING  PERFORMANCE  IS 
EXTREMELY  ROBUST. 


Figure  3.1—1.  Another  unique  aspect  of  nonlinear  control  is  energy  cascade  via  mechanical 
turbulence. 
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impact-damper  control  mechanism  which  received  some  attention  in  the  60’s  to  mid  70’s  (for  exam¬ 
ple,  in  connection  with  buffet  alleviation  in  aircraft,  see  [164,165],  but  which  was  not  subsequently 
pursued  because  of  mechanical  implementation  difficulties.  With  present-day  high-speed  processors, 
however,  such  a  nonlinear  compensator  can  be  implemented  electromechanically  using  a  colocated 
rate  sensor /force  actuator  pair.  However,  it  is  not  suggested  that  research  be  focused  on  the  impact 
idea  per  se.  Rather,  such  devices  are  mentioned  here  solely  to  illustrate  the  potential  of  chaotic 
compensators  and  to  elucidate  some  fundamental  aspects  that  might  be  suitably  generalized  within 
a  rigorous  design  optimization  theory.  The  term  chaotic  compensator  is  used  because,  considered 
by  itself,  the  nonlinear  controller  displays  chaotic  motion.  For  example,  suppose  that  one  discon¬ 
nects  the  chamber  from  the  structure  and  measures  the  response  to  sinusoidal  inputs.  If  there  is 
no  energy  loss  in  collisions,  then  the  system  will  display  homoclinic  tangles  of  great  complexity. 
With  some  energy  loss  mechanisms,  chaotic  attractors  will  result.  Thus,  the  compensator  shown 
in  Figure  3.1-1,  when  considered  alone,  is  a  chaotic  system.  However,  the  intriguing  aspect  here 
is  that  when  this  chaotic  compensator  is  interconnected  with  the  plant,  its  damping  performance 
is  quite  effective  and  extremely  robust.  It  is  important  for  reliable  implementation  of  effective 
compensation  to  understand  and  exploit  the  underlying  mechanisms  involved  in  this  example. 

3.2  Lyapunov  Setting  for  the  Chaotic  Controller 

Lyapunov  theory  provides  the  foundation  for  devising  a  controller  that  emulates  the  behavior 
of  a  chaotic  compensator.  Consider  the  plant  with  dynamics 

*  =  fi(x)  +  /a(*)«» 

V  =  fl  (*)*» 

where  x  £  IRn,  u  €  IRm,  y  £  IRm,  and  dynamic  feedback  controller 

=  /ci  (*« ,  y)  +  As  (*e ,  v)y»  (3) 

«=  -/ea(*«>v)*c  (4) 

where  xe  £  IRn*.  Note  that  the  controller  uses  only  the  available  measurement  y,  although  the 
plant  is  assumed  to  have  a  colocated-type  symmetry  as  in  a  force-to-velocity  model  of  a  flexible 
structure.  We  assume  that  fi(-)  is  dissipative,  that  is, 

xTfi{x)  +  fi{x)x  <  0,  *  €  IRn,  z#0,  (5) 


(1) 

(2) 
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and  that  /i«(*,y)  is  also  dissipative  for  all  y  €  IRm.  In  this  case  the  closed-loop  system  has  the 


£  =  /(*), 


where 


m = 


/l(*)  -  /a(*)/ea(a!e,y)*e 
/« l(*fl.y)  +  /ea(*e.y)/aT(*)* 


with  y  given  by  (2).  Using  the  energy  Lyapunov  function  V(x)  =  zTx  it  is  easy  to  show  that 

-^V(x)  <  0  along  trajectories  of  (6). 
at 

Let  us  now  specialize  to  the  problem  of  vibration  suppression.  Hence  consider  the  plant  model 
i=[-°r2  -2U]*+[i]tl’  * e ntSn,  u e nt1,  (7) 


where 


fl  =  diag  {I?*}  =  modal  frequencies, 
i;  =  diag  {rjk}  =  modal  damping  ratios, 


b  =  modal  actuator  influence  coefficient,  b  €  lRn, 


with  scalar  measurement 


Consider  now  the  compensator 


y  =  bTz3. 


=  -2%]  +2“l°  «Tl*.s)*.  +  *[°]vJ. 


11  =  -k([0  eT]*e)y, 

where  ze  €  IRJn*,  ic  >  0,  o  >  0,  eT  =  [1  1  ...  1], 

fi  =  diag  {/?<},  A  >  0,  i  =  l,.. .  ,nc, 
fj  =  diag  {ft},  ft  >  0,  i  =  l,...,»c, 

and 

■  0  1  1  ...  1- 

-1  0  1  ...  1 

S=  -1  -1  0  ...  1  =_5t. 


-1  -1  -1 
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Note  that  the  plant  and  compensator  are  of  the  form  (l)-(4). 


The  choice  of  feedback  controller  can  be  understood  by  means  of  Figure  3.2-1.  The  underlying 
idea  is  to  transfer  vibrational  energy  from  the  structure  to  the  controller  as  efficiently  as  possible  and 
to  exploit  the  natural  dissipation  of  the  controller.  To  do  this,  the  controller  dynamics  equation  (9) 
involves  an  input  term  proportional  to  y*  to  create  higher-order  harmonics  of  the  natural  structural 
frequencies.  These  harmonics  are  uniformly  distributed  to  a  portion  of  xe  by  means  of  the  vector 
eT  =  [l  1  ...  1].  The  compensator  dynamics  involve  a  dissipative  linear  term 
set  up  its  own  modes  of  vibration.  In  addition,  (9)  involves  a  skew-symmetric  term  S  that  serves 
to  uniformly  distribute,  or  “mix,*  motion  of  all  compensator  states  while  performing  modulation 
(that  is,  creation  of  higher  harmonics)  by  means  of  [0  eT]®c.  Finally,  the  control  signal  given  by 
(10)  again  serves  to  modulate  the  measurement  y  by  the  compensator  harmonics. 


The  intention  of  this  compensator  is  to  purposefully  create  chaos  within  the  controller.  There 
are  two  principal  reasons  for  this  intentional  chaos.  First,  the  structure  itself  has  the  ability  to 
dissipate  energy  by  means  of  the  damping  associated  with  its  natural  modes  of  vibration.  Hence, 
by  creating  higher  frequency  harmonics,  the  compensator  can  efficiently  distribute  low-frequency 
energy,  thereby  exploiting  the  natural  structural  dissipation  to  the  greatest  possible  extent. 


The  second  motivation  for  this  compensator  structure,  as  already  discussed,  is  to  maximize  the 
exchange  of  energy  between  the  plant  and  compensator.  Roughly  speaking,  energy  will  be  trans¬ 
ferred  from  the  structure  to  the  compensator  if  there  is  a  significant  level  of  impedance  matching. 
The  chaotic  motion  within  the  compensator  serves  to  establish  a  broadband  spectrum  to  enhance 
impedance  matching  and  thus  energy  transfer. 

lb  numerically  demonstrate  these  concepts,  we  considered  a  40th-order  (20  modes  between  1 
and  20  rad/sec)  lightly  damped  (.2%  damping)  plant  model  with  a  40th-order  compensator  utilizing 
n  —  n.  To  demonstrate  the  controller  characteristics,  we  considered  the  closed-loop  response  from 
a  nonzero  initial  condition.  Specifically,  the  lowest  frequency  mode  (1  rad/sec)  was  assigned  an 
initial  amplitude  of  unity  and  an  initial  velocity  of  zero,  with  all  other  modes  at  equilibrium. 
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Figure  3.2-1.  The  controller  serves  as  a  mechanism  for  augmenting  energy  dissipation. 
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Figure  3.2-2.  The  time  response  of  the  lowest  frequency  mode  exhibits  rapid  attenuation  due 
to  chaotic  compensator  dynamics. 
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Figure  3.2-2  shows  how  the  amplitude  of  the  first  mode  is  quickly  reduced  to  a  low  level  with  the 
remaining  response  composed  of  broadband  motion.  In  addition,  Figure  3.2-3  shows  the  spectrum 
of  the  measurement  signal  y(t).  This  plot  shows  that  the  structure  undergoes  significant  vibration 
outside  of  the  modal  bandwidth  (approximately  4  Hz).  This  motion,  which  is  due  to  the  nonlinear 
coupling  induced  by  the  controller,  shows  that  energy  is  transferred  from  low  frequency  to  high 
frequency.  Since  the  high  frequency  modes  dissipate  energy  more  efficiently  than  the  low  frequency 
modes  (they  go  to  zero  like  e-'’****),  the  controller  serves  as  an  efficient  mechanism  for  vibration 
suppression. 
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4.0  Phase  Robustness  Theory 

Hoc  theory  accounts  for  modeling  uncertainty  by  bounding  a  weighted  Hoo  norm  characteri¬ 
sation  of  plant  uncertainty.  The  Hoo  norm  does  not  account  for  phase,  however,  which  can  play 
an  important  role  in  robustness  analysis.  For  example,  the  magnitude  of  plant  uncertainty  can  be 
arbitrarily  large  as  long  as  the  phase  of  the  uncertainty  is  such  as  to  avoid  instability. 

Our  interest  in  the  role  of  phase  information  in  robust  control  is  based  upon  connections  to 
power  flow  concepts.  As  will  be  seen,  power  flow  and  stability  analysis  involving  passive  systems  can 
be  extremely  conservative  if  a  small  gain  (Hoo)  approach  is  used.  What  is  lacking  is  the  treatment 
of  phase  properties  which  become  manifested  in  the  structure  of  the  quadratic  Lyapunov  function. 
The  results  described  here  can  be  used  to  guarantee  robust  stability  and  performance  for  both 
linear  and  nonlinear  systems. 

4.1  Positive  Real  Theory  and  Structured  Lyapunov  Functions 

As  a  first  step  in  developing  a  phase  robustness  theory,  we  shall  demonstrate  a  link  between 
phase  properties  and  the  structure  of  the  Lyapunov  function.  Here  we  are  considering  Lyapunov 
functions  of  the  form 

V(x)  =  xTPx  (1) 

where  P  is  a  positive  definite  matrix.  We  shall  call  V(x )  a  structured  Lyapunov  function  if  P  has 
internal  structure.  For  example,  P  may  be  of  the  form 


where  each  diagonal  block  is  also  positive  definite.  We  may,  for  example,  also  require  that  some  of 
the  diagonal  blocks  be  repeated.  Structured  Lyapunov  equations  have  been  studied  in 

S.  Boyd  and  Q.  Yang,  “Structured  and  Simultaneous  Lyapunov  Functions  for  Systems  Stability 
Problems,*  Int.  J.  Contr.  Yol.  49,  pp.  2215-2240,  1990. 

Now  let  us  consider  a  simple  case  of  robustness  due  to  phase.  Consider  the  plant 

x  =  Ax  +  Bu,  (3) 

y  =  Cx,  (4) 
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with  compensator 


xe  =  AcZe  4-  Bcy, 
u  =  —Ccxe. 


(5) 

(6) 


Now  assume  that  the  plant  is  positive  real  and  that  the  compensator  is  strictly  positive  real.  By 
the  Kalman* Yacubovitch  (positive  real)  lemma  there  exist  matrices  L,Le,P,  and  Pc  such  that 


o  =  atp+pa+llt. 

(7) 

PB  =  CT, 

(8) 

0  =  Ajpe  +  PcAe  +  LeLj, 

(9) 

& 

II 

$ 

(10) 

It  is  easy  to  see  intuitively  why  the  closed*loop  system 

~  r-  r  a  f  A  —  BCe 

a=^s<c  A  j 

(11) 

is  asymptotically  stable,  namely,  because  the  phase  shift  of  the  loop  transfer  function  (note  the 
sign  convention  in  (6))  is  less  than  180°.  To  see  this  from  a  Lyapunov  function  perspective,  let  P 

satisfy 

o  =  atp+pa+r, 

(12) 

where 

s  _  Pis 

kTa  * 

(13) 

is  nonnegative  definite.  Expanding  (12)  with 

’  Pi  Pis' 

P  = 

.  ft*  ft  . 

(14) 

yields 

0  =  AtPi  +  PXA  +  (BeC)TP1Ta  -  PxtBeC  +  Ru 

(15) 

0  =  AtPu  +  PnAo  +  (SeC)TPj  -  PiBCc  +  Pia, 

(16) 

0  =  AjP,  +  PaAe  -  {BCe)T P\2  -  P?,BCc  +  P2. 

(17) 

If  we  set 

R\  —  LLT,  Bn  =  0,  R%  =  LeLj 

(18) 
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then  (15)-(17)  are  satisfied  by 


(19) 


Pi  =  P,  Pi J  =  0,  Pj  =  Pc. 


To  see  that  (16)  is  satisfied  note  that  (8)  and  (10)  imply 

( BeC)rPt  -  PiBCe  =  ( BeC)rPe  -  PBCe 


=  CrBjPe  -  PBCe 
=  cTcJ  -  CTCe 


(20) 


=  0. 

Hence  the  Kalman- Yacubovitch  conditions  yield 

which  shows  that  the  Lyapunov  function  for  the  system  “inherits”  the  Lyapunov  function  of  the 
plant  and  compensator. 


P  = 


P  0 
0  Pc 


To  contrast  this  situation  with  Hoo  theory,  suppose  Ru  =  0  but  that  ( BcC)rPt  —  P\BCc  /  0. 
Then  Pit  can  be  suppressed  by  letting 


IIB.II.  ||C.||  «  1 


(22) 


or 

max  Re  A(A«)  <<  0.  (23) 

However,  (22)  and  (23)  correspond  to  small  gain  for  the  feedback  compensator.  The  phase  result, 
however,  does  not  require  either  (22)  or  (23).  Thus  we  have  shown  that  the  Lyapunov  function 
guaranteeing  stability  of  this  feedback  interconnection  has  a  particular  internal  structure.  Since 
the  stability  is  due  to  the  phase  properties  of  the  plant  and  compensator,  we  can  thus  regard  the 
Lyapunov  structure  as  a  manifestation  of  the  phase  aspects. 

4.2  /7-Bound  Theory  and  Structured  Covariances 

Linear  stochastic  control  theory  is  based  on  the  second-moment  statistic  of  the  state  variables. 
Letting  Q  denote  the  state  covariance,  in  the  steady  state  Q  is  given  by  the  Lyapunov  equation 

0  =  AQ  +  QAr  +  V.  (1) 

Suppose  now  that  A  is  uncertain,  that  is,  A  is  replaced  by  A  +  A  A,  where  A  A  €  U,  a  given 
uncertainty  set.  Then  (1)  becomes 

0  =  (A  +  A  A)Qaa  +  Q  ( A  +  A  A)t  +  V.  (2) 
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To  address  (2)  we  introduce  the  notion  of  an  17-bound  which  is  a  matrix  function  satisfying 


A  AQ  +  QAAt  <  17(Q) ,  for  all  Q  >  0,  A  A  ell.  (3) 

That  is,  f](Q)  bounds  the  uncertain  terms  in  (2).  We  now  consider  the  modified  Lyapunov  equation 

0  =  AQ-K3AT+fl(Q)  +  V'.  (4) 

It  is  now  easy  to  show  that  if  (4)  has  a  solution,  then 

Q&.A  <  Q,  for  all  AA  €  U.  (5) 


The  choice  of  17-bound  will  depend  of  course  upon  the  uncertainty  set  U .  However,  for  a  given 
set  U ,  there  may  be  many  17-bounds,  and  a  “best”  bound  need  not  exist  (they  are  only  partially 
ordered).  Two  17-bounds  that  are  convenient  to  work  with  are  the  linear  bound 


and  quadratic  bound 


n{Q)  =  aQ  +  a-l^AiQAj 

»= l 


17(Q)  =  D  +  QEQ. 


(6) 

(7) 


By  choosing  a  special  quadratic  bound,  namely, 


17(Q)  =  ^~3QCTCQ 


(8) 


then  (4)  enforces  an  Hoo  norm  bound  (see  [1-29]). 

The  problem  with  utilizing  bounds  such  as  (6)  or  (7)  is  that  they  may  be  extremely  conservative. 
One  reason  these  bounds  are  conservative  is  that  (3)  must  be  satisfied  for  all  Q  >  0  whether  or 
not  Q  is  the  actual  solution  to  (4).  Moreover,  these  bounds  may  be  conservative  if  the  modeling 
uncertainty  is  large  in  magnitude  but  has  bounded  phase.  Our  approach  to  phase  robustness  theory 
was  motivated  by  the  stochastic  theory  developed  in  [II.1-II.12].  Using  a  multiplicative  noise  model 
with  Stratonovich  interpretation,  Hyland  proposed  the  17-operator 

r 

n{Q)  =  Eli  +  * QA7  +  (Q) 

«=i 

where,  for  a  structural  model  in  modal  coordinates,  each  matrix  A*  is  a  skew-symmetric  matrix 
whose  structure  captures  the  effect  of  an  uncertain  modal  frequency.  A  drawback  of  (9),  however, 
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Quadratic 

Lyapunov  Functions 


and  unstructured,  complex  uncertainty. 


Structured  Lyapunov 
Function 


Figure  4.2  2.  Phase  information  is  manifested  in  structured  Lyapunov  functions  and  is  directly 
related  to  structure,  real-valued  uncertainty. 


is  that  fi(Q)  is  indefinite.  Thus,  in  this  case  the  modified  Lyapunov  equation  (4)  does  not  provide 
a  bound  for  Qaa  and  thus  does  not  guarantee  stability  by  means  of  standard  techniques. 

In  summary,  we  note  that  there  is  an  intricate  interplay  between  phase  information,  real  par 
rameter  uncertainty,  and  Lyapunov  functions.  The  classical  situation  shown  in  Figure  4.2-1  is 
thus  an  extreme  case  of  the  more  subtle  situation  addressed  in  Figure  4.2-2.  Further  discussion  of 
these  issues  can  be  found  in  the  paper  “Real  Parameter  Uncertainty  and  Phase  Information  in  the 
Robust  Control  of  Flexible  Structures,”  which  appears  in  Appendix  E. 

4.3  17-Bounds  for  Positive  Real  Theory 

To  exploit  the  features  of  positive  real  transfer  functions,  we  have  developed  a  theory  of  robust 
controller  synthesis  with  positive  real  uncertainty.  The  phase-bounded  character  of  positive  real 
transfer  functions  entails  far  less  conservatism  than  small  gain  or  H,*  results  when  addressing  real 
parameter  uncertainty. 

The  results  obtained  thus  far  are  detailed  in  the  paper  “Robust  Stabilization  with  Positive 
Real  Uncertainty:  Beyond  the  Small  Gain  Theorem,”  which  appears  in  Appendix  F.  This  paper 
develops  a  state  space  theoiy  of  positive  real  transfer  functions  in  terms  of  an  algebraic  Riccati 
equation.  This  characterization  is  more  direct  than  the  usual  KYP  characterization  and  provides 
the  basis  for  state  space  controller  synthesis  techniques  in  the  spirit  of  state  space  Hqo  theory. 

More  recently  we  have  linked  positive  real  theory  with  /7-bound  theory  by  showing  that  robust 
stability  and  robust  Hoo  performance  in  the  presence,  of  positive  real  uncertainty  are  guaranteed  by 
means  of  an  /7-bound.  This  connection  has  ramifications  for  nonlinear  control.  To  see  this,  we  re¬ 
call  that  robust  stability  in  the  presence  of  sector-bounded  nonlinearities  is  equivalent  to  a  Nyquist 
circle  criterion,  which  is  equivalent  to  a  positive  real  condition.  Thus  robustness  to  positive  real  un¬ 
certainty  provides  the  means  to  guarantee  stability  with  respect  to  a  class  of  nonlinearities.  Similar 
observations  hold  for  the  Popov  criterion  which  also  guarantees  robustness  for  sector  nonlinearities. 

Our  results  provide  the  means  for  going  beyond  existing  results  in  two  respects.  First,  we  can 
develop  multivariable  generalizations  of  the  classical  circle  and  Popov  criteria  using  simplified  /7- 
bound  theory.  And,  second,  our  techniques  can  be  used  for  robust  synthesis  in  addition  to  analysis 
as  addressed  by  standard  theory. 
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5.0  Optimal  Nonlinear  Feedback  Control 


The  methods  and  results  discussed  in  Sections  2-5  are  independent  of  optimality  considerations. 
The  purpose  of  this  section  is  to  discuss  progress  in  developing  an  optimality-based  control  theory 
involving  nonlinear  controllers  for  linear  and  nonlinear  plants. 

As  pointed  out  in  Section  2,  controller  synthesis  need  not  be  based  upon  optimality  crite¬ 
ria.  For  example,  a  controller  can  be  constructed  in  accordance  with  a  Lyapunov  function  to 
achieve  stability,  energy  dissipation,  etc.  Nevertheless,  there  is  strong  motivation  for  developing  an 
optimality-based  theory. 

Perhaps  the  prime  motivations  for  developing  an  optimality-based  theory  is  the  success  of 
such  approaches  in  linear  control  theory.  The  well-known  linear-quadratic-Gaussian  control  theory 
(LQG)  is  the  major  result  in  modem  optimal  multivariable  feedback  control  theory.  During  the 
past  decade,  LQG  theory  has  been  extended  to  address  numerous  practical  control-design  issues 
such  as  disturbance  attenuation,  robustness,  controller  order,  and  pole  placement  (Figure  5-1).  The 
resulting  theory,  known  as  Optimal  Projection  for  Uncertain  Systems  (OPUS),  has  been  extensively 
developed  (see  the  reference  list  in  Appendix  B). 

The  second  motivation  for  optimal  nonlinear  control  theory  is  that  it  can  drive  the  controller 
synthesis  procedure  within  a  class  of  candidate  controllers.  Specifically,  as  will  be  discussed  later 
in  this  section,  we  can  view  a  given  Lyapunov  function  as  providing  the  framework  for  controller 
synthesis  by  guaranteeing  local  or  global  asymptotic  stability  theory  for  a  class  of  feedback  con¬ 
trollers.  The  actual  controller  choeen  for  implementation  can  thus  be  the  member  of  this  candidate 
class  that  minimises  a  specified  performance  function.  The  form  of  this  functional  is  usually  closely 
related  to  the  structure  of  the  Lyapunov  function.  In  LQG  theory,  for  example,  the  Lyapunov 
function  is  the  familiar  quadratic  function  V(z)  =  xTPx,  while  the  gains  are  chosen  to  minimise 
a  performance  functional  of  the  form  J  =  tr  PV .  In  summary,  then,  Lyapunov  function  theory 
provides  the  framework ,  while  optimisation  fixes  the  gains. 

5.1  Optimal  Nonlinear  Feedback  Control  via  Steady-State  HJB  Theory 

The  classical  approach  to  optimal  nonlinear  control  is  to  invoke  the  Maximum  Principle.  This 
result  has  been  successful  in  characterizing  solutions  to  problems  such  as  minimum  time  con¬ 
trol.  Since  the  Maximum  Principle  does  not  explicitly  guarantee  stability  via  a  Lyapunov  function 
per  se  and  does  not  directly  lead  to  feedback  controllers,  we  shall  not  adopt  it  as  our  principal 
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Figure  5—1.  Optimal  Projection  for  Uncertain  Systems  (OPUS)  is  an  optimal  linear  control  theory 
that  systematically  addresses  a  broad  range  of  practical  control  design  issues. 
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approach.  Nevertheless,  we  bear  in  mind  that  the  Maximum  Principle  does  have  links  with  the 
Hamilton- Jacobi- Be  liman  (HJB)  approach  which  we  shall  consider  and  to  which  we  now  turn. 

Hamilton- Jacobi- Bellman  theory  has  its  roots  in  the  classical  Hamilton- J acobi  partial  differen¬ 
tial  equation  as  well  as  the  dynamic  programming  technique  of  Bellman.  In  its  most  general  form, 
the  theory  involves  a  partial  differential  equation  whose  solution  yields  an  optimal  controller.  In 
recent  years,  this  equation  has  attracted  renewed  interest  with  the  discovery  of  generalised  solutions 
[151,152]. 

If,  in  accordance  with  practical  motivations,  we  restrict  our  attention  to  time-invariant  systems 
on  the  infinite  horizon  with  analytic  data,  the  situation  is  considerably  simplified.  In  this  case  the 
HJB  partial  differential  equation  reduces  to  a  purely  algebraic  relationship. 

lb  summarise  the  ideas  involved  we  first  consider  the  problem  of  evaluating  a  nonquadratic  cost 
functional  depending  upon  a  nonlinear  differential  equation.  It  turns  out  that  the  cost  functional 
can  be  evaluated  in  closed  form  so  long  as  the  cost  functional  is  related  in  a  specific  way  to  an 
underlying  Lyapunov  function.  The  basis  far  the  following  development  is  the  paper  [60]  by  Bass 
and  Weber.  A  more  detailed  treatment  of  these  results  is  given  in  the  paper  “Nonquadratic  Cost 
and  Nonlinear  Feedback  Control”  which  appears  in  Appendix  G. 

For  simplicity  in  the  exposition  here,  we  shall  define  all  functions  globally  and  assume  that 
existence  and  uniqueness  properties  of  the  given  differential  equations  are  satisfied. 

For  the  following  result,  let  /:  IR"  — ►  IR"  and  L:  IR"  — *  IR.  We  assume  /(0)  =  0. 

Lemma  1.  Consider  the  system 

*(<)  =  /(*(*))>  *(°)  =  *o,  (1) 

with  performance  functional 

•/(*,)=  (2) 

Jo 

Assume  that 

L(x)  >  0,  x  €  IR",  x  #  0,  (3) 

and  assume  there  exists  a  Cx  function  V:  IR"  — »  IR  such  that 

V(0)  =  0,  (4) 

V(x)>0,  x  €  IR",  x#0,  (5) 

L(x)  =  -H*)/(*),  *€IR".  (6) 


H»rr is  Corp. 


5-3 


Daetnbu  1090 


Then  x  =  0  is  a  globally  asymptotically  stable  solution  of  (1)  and,  furthermore, 


J(x0)  =  V(x0). 


(7) 


Proof.  Let  x(t)  satisfy  (1).  Then 


v(«M)  =  |v(*(.))  = 


(8) 


Hence  it  follows  from  (6)  that 


Now  (3)  implies  that 


V(*(t))  =  -L(*(t)). 


V(x(t))<  0,  x(t)  ?  0. 


Since  V(x)  >  0,  x  #  0,  it  follows  that  V(-)  is  a  Lyapunov  function  for  (1)  and  that  x(t)  — ♦  0  as 
t  — ►  oo.  Thus  proves  global  asymptotic  stability  of  the  solution  x  =  0.  Now  (8)  implies  that 

v'WO)  -  v(*>)  =  /*v'(x(.))/(x(.)) d. 

Jo 

=  -  /  £(*(«))ds. 

Jo 

Letting  t  —*  oo  and  noting  V(x(t))  — >  0,  it  follows  that 

-V (xo)  =  -  r  L(x(t))dt, 

Jo 

or,  equivalently, 

V(x0)  =  J(x0).  □ 


The  main  feature  of  Lemma  1  is  the  role  played  by  the  Lyapunov  function  V  (x)  in  guaranteeing 
stability  and  for  evaluating  the  functional  J(x o).  It  can  be  recognised  that  V(x)  is  the  cost- to- go 
function  in  dynamic  programming. 

Let  us  illustrate  Lemma  1  with  a  familiar  example.  Consider  the  linear  system 

x  =  Ax,  z(0)  =  x0,  (9) 

with  cost  functional 

J(x o)  =  f  xT Rx  dt,  (10) 

Jo 
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where  R  €  IR"X"  is  positive-definite.  If  A  is  stable  then  there  exists  a  positive  definite  matrix 
P  €  TRnxn  satisfying 

0  =  ArP  +  PA+R.  (11) 


Now  define 

V(x)  =  xTPx,  (12) 

which  satisfies  (4)  and  (5).  Furthermore,  with  /(x)  =  Ax  and  L(x)  =  xT Rx  it  follows  that 

V#(x)/(x)  =  2  xrPAx 

=  xT(ATP  +  PA)x 
=  -xTRx 
=  -Hx) 

which  verifies  (6).  Hence 

J(x0)  =  XqPxq, 

which  is  a  familiar  result  from  linear-quadratic  theory. 


To  deal  with  more  general  situations,  the  following  lemma,  which  appears  in  [60],  will  be  useful. 

Lemma  2.  Let  A  €  lRnXn  be  asymptotically  stable  and  let  h:  1R”  -+•  IR  be  a  nonnegative- 
definite  homogeneous  p-form  (p  even).  Then  there  exists  a  nonnegative-definite  homogeneous  p-form 
g:  IRn  — » 1R  such  that 

g'(x)Ax  +  h(x)  =0,  x  G  St".  (13) 


In  the  quadratic  case  (13)  yields  the  familiar  result.  To  see  this  let  h{x)  =  xTRx  and  g(x)  — 
xTPx.  Then  (13)  becomes 

2xTPAx  +  xtRx  =  0, 


or 


xT(ATP  +  PA  +  R)x  =  0, 

which  is  satisfied  by  P  given  by  (11).  Now  consider  the  nonquadratic  cost  functional 

J(xo)  =  [  [xtRx  -f  h(x)]dt, 

Jo 

where 

M*)  =  5Z 


(14) 


(15) 


v-\ 
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,  r,  /ij„:  IRn  — ► IR  is  a  nonnegative-defmite  homogeneous  2i/-form.  We  continue 


and,  for  v  =  1, . . . 

to  assume  that  x(t)  satisfies  (9),  where  A  is  stable.  Now,  let  gjv:  lRn  — *  IR  be  the  nonnegative- 
definite  homogeneous  2t/-form  satisfying 

g3u{x)Ax  +  ha„(x)  =  0,  x  e  IRn,  v  =  1, . . . , r,  (16) 

and  define 

r 

«(*)  =  £«  IvOO-  (17) 

K=1 

Note  that  (15)-(17)  imply 

g'(x)Ax  +  h(x)  =  0,  x  €  IRre.  (18) 

Furthermore,  define  the  positive-definite  function 

V(x)  =  xTPx  +  g(x),  (19) 

where  P  satisfies  (11).  Now,  to  verify  (6)  we  note  that 

V'(x)/(x)  =  (2xtJP  +  y'(x)]Ax 

=  XT(ATP  +  PA)x  +  ^2  g'3v(x)Ax 
v-\ 

=  -xtRx  -  ^2  W*) 

y=l 

=  -L(x). 

Hence  for  J(x0)  given  by  (14)  we  obtain 

J  (*o)  =  V'(xo)  =  xj  Px  o  +  g{x0) .  (20) 

Next  consider  in  place  of  (9)  the  case  in  which  the  plant  is  nonlinear,  for  example, 

x  =  Ax  +  <r{x),  x(0)  =  xo,  (21) 

where  <r(0)  =  0  and  we  continue  to  assume  that  A  is  stable.  Again,  let  g:  lRn  — ►  IR  be  given  by 
(16)  and  (17)  and  define  V(x)  by  means  of  (19).  It  remains  only  to  verify  (6).  Hence 

v'{x)f(x)  -  [2 xTP  +  j'(x)][Ax  +  a(x)] 

=  xT(ATP  +  PA)x  +  g,(x)Ax  +  [2  xTP  +  g'{x)\a{x) 

=  -[xTPx  +  h(x)]  +  [2xTP  +  p#(x)]a(x) 

=  -{£(*)  “  (2xTP  +  $'(x)]<r(x)}. 

5-6 
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Hence  we  see  that  (6)  is  not  generally  satisfied.  However,  we  can  salvage  the  situation  be  considering 
an  auxiliary  cost  functional 

i(*o)=  rt(x(i))it,  (22) 

Jo 

where 

L(x)  4  L(x)  -  [2 xTP  +  g'{x))a{x).  (23) 

With  this  modification  (3)  must  be  satisfied  with  L(x)  replaced  by  L(x),  that  is, 

L(x)  >  [2xT  P  +  g'(x)]a(x)  (24) 

In  the  special  case  that 

[2*TP  + j'(*)]o(i)  <  0,  (25) 

it  follows  that  (24)  is  automatically  satisfied  (if  (3)  is  satisfied)  and,  furthermore, 

J(x0)  <  J(*0).  (26) 

That  is,  the  auxiliary  cost  is  an  upper  bound  for  the  original  cost. 

By  only  a  slight  extension  of  Lemma  1,  we  can  obtain  sufficient  conditions  for  characterizing 
optimal  feedback  controllers.  Now  let/:  IRnxIRm  — »  IR",  where/ (0,0)  =  0,  let  L:  IRnxlRm  — ♦  IR, 


and  define  for  p  €  IRn 

H{x,p,u)  =  L(*,u)  +pT/(x,u). 

Theorem  1.  Consider  the  controlled  system 

*(0  =  /(*(*)>«(*))•  *(0)  =  x0,  (27) 

with  performance  functional 

J(x0, «(•))-  f  L(x(t),u(t))dt.  (28) 

Jo 

Assume  that  there  exist  a  Cl  function  V :  IR"  — *  IR  and  a  function  <f>:  IR"  — ♦  IRm  such  that 

V(0)  =  0  (29) 

V(x)  >  0,  x  G  IR",  x^0,  (30) 

L(x,^(x))>  0,  xeIR",  x#  0,  (31) 

H{xtVn(x)A{ *))  =  0,  x  €  IRn,  (32) 

H(x,V'r{x),u)  >0,  z€lR",  u  €  IRm.  (33) 
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Then  with  the  feedback  control  u(-)  =  ^(z(>)),  the  solution  z  =  0  of  the  closed-loop  system  is 
asymptotically  stable  and 

J(x0, *(*(■))  =V(x0).  (33) 

Furthermore,  the  feedback  control  u(-)  =  <fr(x(-))  minimizes  J(x o,u(*)),  that  is, 

J(*o,^(*(-))  =  min  J(*0,  «(•))•  (34) 

«(•) 

Proof.  Global  asymptotic  stability  and  result  (33)  follow  directly  from  Lemma  1.  We  need 
only  note  that  (31)  can  be  written  as 

L(x,  *(x))  =  -V(x)/(x,  4(x)),  x  €  IRn, 

which  corresponds  to  (6).  It  remains  only  to  prove  (34)  using  condition  (32).  For  arbitrary  u(t) 
and  for  x(t)  satisfying  (15)  we  have 

V(x{t))=V'{x(t))f(x{t)Mt)) 

or 

0  =  -V(x(t))  +  V'(x(t))f{x{t),u(t)). 

Hence 

£(*(*). «(0)  =  +  £(*(<).«(<) + VW0/W0.-W) 

=  -V(x(()  +  *(*«,  v'T(x(f)),«(t)). 

Now  using  (32)  and  (33)  we  obtain 

J(*o,«(0)  = 

Jo 

=  -  lim  K(z(t))+V(z0)  +  /°°H(*(t),^/T(*(0),«(0)d« 

*-*oo  J0 

-V(x0)  +  [  *(*(*)»  ^**(*(0*  «(*))** 

Jo 

>  V{xo) 

=  J[x  0,^(x(-)), 

which  yields  (33).  □ 

The  principal  feature  of  Theorem  1  is  that  the  optimal  control  law  u  =  4>(x)  is  a  feedback 
controller.  Furthermore,  this  control  is  optimal  independently  of  the  initial  condition  zq. 
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Now  let  us  illustrate  Theorem  1  with  some  examples.  We  begin  with  the  simplest  case,  namely, 
the  linear  quadratic  regulator.  Hence  consider  the  controlled  system 

x  =  Ax  +  Bu,  x(0)  =  *o.  (36) 

with  performance  functional 

J(x0,  «(•))=  f  [xTi2ix  +  uTi?ju]dt,  (37) 

Jo 

where  x  €  IRn,u  €  IRm  and  where  R\  and  Ri  are  positive  definite.  Define  the  feedback  law 

*(x)  =  -RjlBTPx,  (38) 

where  P  satisfies 

0  =  ArP  +  PA +  Hi  -  PSP ,  (39) 

where  S  =  BRflBT.  With  u  =  ^(x),  the  closed-loop  system  (36)  becomes 


x  =  Ax,  x(0)  =  xq, 


(40) 


(41) 

(42) 


where  R  =  Ri  +  PSP.  Thus  the  closed-loop  system  (40)  with  cost  (41)  has  exactly  the  form  of 
the  example  considered  in  (9)-(12).  It  remains  only  to  show  that  u  =  <f>(x)  is  the  optimal  control, 
which  will  be  the  case  if  (33)  is  satisfied.  To  show  this,  note  that 


£T(x,  V,T(x),u)  =  xrRxx  +  uTfZju  +  2xrP(Ax  -I-  Bu) 

—  xTI?ix  +  uTiiju  +  xt(AtP  +  PA)x  +  2  xTPBu 

—  xTPSPx  +  2xrPBu  -(-  uTi?ju 

=  [R^xBtPx  +  u)T  R2[Rjl  Br  Px  +  u] 

>  0. 


Note  that  it  is  now  easy  confirm  (32)  by  setting  u  =  4>{x)  =  -Hf 1  Br Px. 
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We  now  apply  Theorem  1  to  an  optimal  control  problem  involving  a  nonquadratic  cost.  Hence 
consider  the  linear  system  (35)  with  cost 

J (*o, «(•))  =  [  [xtR\X  +  h(x)  +  uT.Rju]dt,  (43) 

Jo 

where  h(x)  is  given  by  (15).  We  shall  consider  a  control  law  of  the  form 

u  =  4(x)  =  4>it(x)  +  <f>Ni.{x)>  (44) 

where  ^r,(x)  and  4nl{*)  are  linear  and  nonlinear,  respectively.  Let  4>l{x)  agree  with  the  linear- 
quadratic  solution,  that  is, 

<f>L{x)  ±  -R^lBrPx,  (45) 


where  P  satisfies 


0  =  ArP  +  PA  +  Ri-  PSP. 


Recall  that  (46)  can  be  written  as 


0  =  AtP  +  PA+R, 


where  A  =  A-  SP  and  R  =  Ri  +  PSP. 


For  the  nonlinear  control  ^jvl(*)  let  $(x)  be  given  by  (17)  where  03„(z)  satisfies 

g'iv(x)Ax  +  h7v(x),  x  e  IRn,  v  =  1 . . 

which  is  the  same  as  (16)  with  A  replaced  by  A.  Now  define 

4>nl{?)  =  -!flf1£VT(z) 


and  the  Lyapunov  function 


V(x)  =  xTPx  +  g(x). 


Next  note  that  for  L(x,  u)  as  in  (43)  we  have 


L(x,f{x))  =  xrRix  +  h(x)  +  ^T(x)f?2^(z) 

=  xtRx  +  h(x)  +  xTPSg,T(x)  +  lg'(x)Sg'T(x). 

4 


Furthermore,  with  u  =  (f>{x)  the  system  (36)  becomes 


x  =  Ax  +  B<I>sl{x ),  x(0)  =  xo, 
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or 

x  =  Ax-  ^Sg'T(x),  i(0)  =  x0.  (53) 

Returning  to  Theorem  1,  it  is  clear  that  (29)-(31)  are  satisfied.  However,  note  that 
V '(*)/(**  ^(*))  -  (2 xrP  +  g'(z)][Ax  -  -Sff/T(x)j 

=  xr(ATP  +  PA)x  +  y'(x)Ax  -  xTPSg'T(x)  -  \g'(x)Sgn:(x) 

It 

-  -  [xtRx  +  h(x)  +  xTPSg,T(x)  +  ^g'(x)5g,T(x)] 

=  -  [£(**#*))  +  ^s'(*)VT(*)]> 

so  that 

=  -jff'WSg^Cz),  (54) 

which  shows  that  (32)  is  not  satisfied.  However,  if  we  define 

L(x,  u)  =  L(x,  u)  +  ^g'(x)SgIT(x),  (55) 

then  the  auxiliary  cost 

•J(*o,  «(•))-  f  L(x(t),u(t))dt  (56) 

Jo 

satisfies 

J(x0,u(-))<i(xo,  «(•)).  (57) 

Finally,  be  defining 

H(x,Vn  (x),u)  =  L(x,  u)  +  Vn(x)f(x,u),  (58) 

it  can  be  shown  that 

H(xtVn(x)tu)  =  [u  -  ^(x)]TRa[u  -  4>(x)\.  (59) 

Hence  (33)  holds  with  Zf(-)  replaced  by  #(■).  Consequently, 

7(x0,^(x(-))  =  V(x0) 

=  xj  Px  o  +  g(x0)  (60) 

=  nun  7(x0,  «(•))• 

u() 

We  next  consider  a  special  case  of  the  above  nonquadratic  problem  that  leads  to  considerable 
simplification.  This  particular  problem  was  considered  in  [63].  Suppose  we  require  that  V(x)  be  of 

the  form 

V(x)  =  xTPx  +  i(xTMx)3  (61) 

mt 
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so  that  g(x)  =  ^(iTAfi)1  where  P  satisfies  (39)  and  M  satisfies 


0=  (A  -  SP)tM  +  M(A  -  SP )  +  #!  +  MSM.  (62) 


Then  4(x)  has  the  form 


4>{x)  =  -RflBTPx  -  R}1  BT  (xT  Mx)Mx. 


(63) 


Next  we  assume  that  h(x)  is  given  by 

h(x)  =  (xTMx)xT(Ri  +  MSM)x.  (64) 

With  these  definitions  we  note  that 

g'(x)Ax  =  (xT  Mx)2xr  MAx 

=  (xtMx)xt(AtM  +  MA)x 
=  (xtMx)xt(R1  +  MSM)x 

=  M*)> 

which  verifies  (48).  Finally,  define 


L[x,  u)  =  x"^  R\X  +  h{x).  (65) 

Following  the  previous  development,  we  see  that  ^(z)  given  by  (62)  minimizes  J(x0,  u(-))  defined 
by  (55),  where 

L(x,  u)  =  xtR\X  +  (zTilfz)zT(fJi  +  MSM)x  +  (xT  Mx)*xr  MS  Mx  +  uTR2u.  (66) 

Thus,  be  minimizing  a  sixth-order  cost  functional,  the  optimal  control  is  a  cubic  feedback  charac¬ 
terized  by  a  pair  of  Riccati  equations.  The  cost  functional  is  somewhat  artificial  since  it  depends 
upon  the  solution  of  one  of  the  Riccati  equations. 

We  have  thus  shown  that  the  problem  considered  by  Speyer  in  [63]  is  a  special  case  of  the 
optimal  nonquadratic  cost  problem  addressed  by  Bass  and  Weber  in  [60].  Actually,  however,  the 
formulation  of  Speyer  was  a  stochastic  control  problem  based  upon  results  of  Wonham  [153].  This 
formulation  involves  systems  of  the  form. 

x  =  Ax  +  Bu  +  Diw,  (67) 


Harris  Corp. 


5-12 


December  1990 


when  D\w  denotes  additive  white  noise  disturbances.  (Speyer  also  considered  multiplicative  noise 
in  [63]  as  well.)  These  disturbances  lead  to  a  modification  of  (46)  of  the  form 

0  =  ArP  +  PA  +  Rt-  PSP  +  (tr  MVt)M  +  2MViM,  (68) 

Now  note  that  (62)  and  (68)  now  constitute  a  pair  of  coupled  Riccati  equations. 

Having  reviewed  the  elements  of  a  deterministic  HJB  theory  as  originated  by  Bass  and  Weber, 
out  next  goal  is  to  develop  a  corresponding  theory  of  stochastic  control.  Such  a  theory  can  be  used 
for  disturbance  rejection  for  persistent  disturbances.  Our  principal  goal,  however,  is  to  generalize 
HJB  theory  to  permit  the  design  of  fixed-structure  controllers  that  operate  on  the  available,  possibly 
noisy,  measurements.  To  our  knowledge,  no  such  theory  currently  exists,  while  progress  in  this 
direction  is  crucial  for  practical  application  of  nonlinear  control  laws. 
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Conf.,  pp.  1590-1597,  Seattle,  WA,  June  1986. 
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H.38  D.  S.  Bernstein,  “Optimal  Projection/Guaranteed  Cost  Control  Design  Synthesis:  Robust 
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11.47  F.  M.  Ham  and  S.  W.  Greeley,  “Active  Damping  Control  Design  for  the  Mast  Flight  System,” 
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Net.  Sys.,  Phoenix,  AZ,  June  1987. 
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August  1987. 
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H.65  D.  C.  Hyland,  “Experimental  Investigations  in  Active  V'Sration  Control  for  Application  to 
Large  Space  Systems,”  in  Space  Structures,  Power,  and  Power  Conditioning,  R.  F.  Askew, 
Ed.,  Proc.  SPIE,  Vol.  871,  pp.  242-253,  1988. 
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Abstract 

It  is  well  known  from  thermodynamics  that  energy  flows  from  hot  objects  to  cold 
objects.  It  is  less  well  known,  however,  that  a  similar  phenomenon  occurs  in  coupled 
mechanical  systems  with  modal  energy  playing  the  role  of  temperature.  Energy  flow  among 
coupled  modes  is  the  subject  of  Statistical  Energy  Analysis  (SEA).  Originally  motivated  by 
problems  in  acoustics  involving  numerous  vibrational  modes,  SEA  is  based  upon  equations 
governing  energy  flow  among  individual  modes  or  sets  of  modes.  Such  energy  flow  equations 
can  be  quite  efficient  in  modeling  the  response  of  lightly  damped  structures.  This  paper 
has  two  goals.  First,  we  derive  a  generalized  formulation  of  power  flow  which  allows 
arbitrary  coupling  of  arbitrary  strength.  Previous  theoretical  results  were  limited  to  either 
identical  couplings  or  weak  interactions.  These  new  results  utilize  Kronecker  matrix  algebra 
to  derive  an  energy  flow  equation  involving  the  diagonal  elements  of  the  solution  to  a 
Lyapunov  equation.  Analysis  of  the  resulting  equations,  based  upon  M-matrix  theory, 
yields  generalized  energy  balance  relations  in  the  case  of  weak  but  arbitrary  (possibly 
nonconservative)  couplings. 
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Notation 


IE 

IR,<D 

lRrx*,  (Drx* 
IRr,Cr 
Ir  or  I 
3 

Ant- 

Re  A,  Im  A 
A,  AT,  A* 

{A},  (A) 

®,©,  vec.vecd 
A  >>  0 


expectation 

real  field,  complex  field 
r  x  a  real,  complex  matrices 
IRrXl,Crxl  (column  vectors) 
r  x  r  identity  matrix 

V=T 

(Jfc,  ^-element  of  A  6  Crx* 
real,  imaginary  part  of  A  €  <D,X* 

complex  corrugate,  transpose,  complex 
conjugate  transpose  of  A  €  <Crx* 

diagonal,  off-diagonal  part  of  A  €  Crxr  (see  Section  2) 

See  Appendix  B 

A  €  IRrX*  is  nonnegative  (each  entry  of  A  is  nonnegative) 
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1.  Introduction 


We  are  concerned  with  efficient  methods  for  evaluating  the  steady-state  statistical  response  of 
large-scale  linear  systems  composed  of  many  interconnected,  high-dimensioned  subsystems.  This 
problem  arises  from  applications  involving  acoustical  response,  acoustical  /structural  interaction, 
high  frequency  vibration  of  mechanical  systems,  and  dynamics  and  control  of  large  space  struc¬ 
tures  [1-25].  To  illustrate  the  problem,  suppose  that  each  subsystem  is  well  known  and  precisely 
characterized  so  that  its  eigenbasis  is  known.  Then,  assuming  for  convenience  a  semisimple  eigen- 
structure,  the  A:th  subsystem  considered  in  isolation  is  of  the  form 

ik  =  +  wk,  k=  l,...,r,  (1.1) 

where,  for  k  =  1 . . 

*fc€Cn\  Afe=  diag  (A*),  A*  €  <D, 

♦=i, —.*»* 

and  wjb  is  a  white  noise  process  with  Hermitian  nonnegative-definite  intensity  Vk  e  <DnfcXnfc.  When 
the  subsystems  are  interconnected,  couplings  are  introduced  among  the  subsystems  in  the  form  of 
perturbations  to  the  individual  subsystems.  The  subsystem  dynamics  in  the  interconnected  case 
are  given  by 

r 

xk  =  Afcgfc  -+  gkkxk  +  yi 9Uxt  +  wk>  k=  1,. . . ,r.  (1.2) 

lm  1 

The  matrix  gkt  e  <D"*Xn<,  k  ^  l,  represents  the  effect  of  the  Zth  subsystem  on  xk,  while  the 
matrix  gkk  €  <Cnfc  Xn*  represents  an  effective  shift  of  A*  due  to  the  interconnections.  Our  goal  is  to 
determine  the  steady-state  second-moment  response  of  the  interconnected  systems. 

In  the  case  of  a  large  flexible  structure  consisting  of  several  interconnected  substructures, 
(1.1)  represents  the  Ath  substructure,  while  the  coupling  terms  in  (1.2)  arise  from  the  mechani¬ 
cal  interconnections  among  the  substructures.  Alternatively,  in  the  case  of  acoustical/structural 
interaction,  one  wishes  to  predict  the  acoustical  response  of  several  acoustical  spaces  separated 
by  elastic  partitions  (such  as  walls).  Equation  (1.1)  then  represents  the  modal  dynamics  of  each 
acoustical  chamber,  while  the  coupling  terms  in  (1.2)  represent  the  dynamic  couplings  introduced 
by  the  Jastic  partitions.  The  power  flow  concept  also  has  close  connections  with  thermodynamics 
and  circuit  theory  [26-34] . 

The  problem  posed  by  (1.2)  is  subsumed  in  the  linear  system  model 

x  —  {-v  +jn)x  -t-  (H  +  G)x  +  w,  (1.3) 
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where 

x  e  cn, 

v=  diag  (i/fc),  I/Jb  €  IR,  Vk  >  0,  k=l,...,n, 

k=l,...,n 

n  =  diag  (/?*),  Hk  €  IR,  fc  =  l,...,n, 

. » 

diag  (fffc),  #*€<0,  fc=l,...,n, 

*=i . *» 

Gfcfc  =  0,  *  =  l,...,n,  G  €  <DnXn, 

and  w  is  white  noise  with  Hermitian  nonnegative-definite  intensity  V  6  CnXn.  The  diagonal 
matrix  —  v  +  jf2  €  <Cnxn  where  n  =  St=i  nt  is  a  concatenation  of  all  of  the  uncoupled  subsystems 
in  (1.1).  The  matrices  H  and  G  represent,  respectively,  the  diagonal  and  off-diagonal  portions  of 
the  perturbations  due  to  subsystem  interaction.  We  asstime  that  the  system  (1.3)  is  asymptotically 
stable,  that  is,  the  spectrum  of  the  matrix  —  +  H  +  G  ia  contained  in  the  open  left  half  plane. 

To  study  the  steady-state,  mean-square  response  of  the  system  (1.3),  suppose  y  defined  by 

y  —  Cx  (1.4) 

is  a  response  variable  of  interest,  where  C  =  [C\  ••• Cn\  €  <Dlxn.  Then  it  is  well  known  [35]  that 
the  steady-state  mean-square  value  of  y  is  given  by 

limE[|y]3]  =  tr[C*CQ],  (1.5) 

where  the  steady-state  covariance  Q  =  limt_oo  IEfarrc*]  6  <CnXn  is  determined  as  the  Hermitian 
nonnegative-definite  solution  to  the  Lyapunov  equation 

0  =  (-1/  +jn)Q  +  Q{-u  -jn)  +  {H  +  G)Q  +  Q(H  +  G)*  +  V.  (1.6) 

Note  that  due  to  symmetry,  equation  (1.6)  represents  -§n(n  +  1)  scalar  equations  for  the  elements 
of  Q. 

Since  (1.6)  is  a  well-known  equation  with  well-established  solution  techniques  [36-38],  the 
problem  would  appear  to  be  solved.  However,  the  difficulty  in  the  application  mentioned  above 
is  that  the  total  system  dimension  n  may  be  exceedingly  large.  For  the  example  involving  several 
acoustic  spaces  coupled  by  elastic  partitions,  each  subsystem  (a  modest-sized  room,  say)  can  have 
millions  of  modes  in  the  audio  range.  Thus  the  total  dimension  n  can  be  of  the  order  of  10e  —  107 
while  the  coefficient  matrix  (-i /  +jfl  +  H  +  G)  is  not  necessarily  either  sparse  or  banded.  Thus, 
the  prediction  of  vibration  response  or  sound  pressure  levels  via  the  solution  of  (1.6)  can  be  very 
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cumbersome  indeed.  It  is  thus  desirable  to  develop  more  efficient  methods  for  estimating  quantities 
such  as  IE[|y|3]  which  somehow  circumvent  the  huge  dimensionality  of  (1.6). 

In  this  regard,  many  useful  and  important  results  and  procedures  have  been  developed.  These 
are  often  referred  to  collectively  as  “Statistical  Energy  Analysis,”  or  SEA  [1-16].  SEA  was  developed 
for  high-dimensional,  lightly  damped  mechanical  or  acoustical  systems  for  which  there  are  passive 
mechanical  energy-conservative  interconnections  among  the  subsystems.  In  the  notation  of  (1.3), 
(1.4),  this  means  that  there  is  a  basis  in  which  H  -f  G  is  skew-Hermitian,  that  is, 

Hk=jSk ,  Hk  €  IR,  &=l,...,n,  (1.7) 

G*  =  —G.  (1.8) 

In  the  present  paper,  we  develop  results  that  deal  with  general  coupling  terms.  These  results  are 
later  specialized  to  couplings  restricted  by  (1.7),  (1.8). 

The  purpose  of  this  paper  is  to  elucidate  some  of  the  basic  ideas  of  SEA  in  rigorous  system- 
theoretic  language  and  to  provide  generalizations  of  certain  fundamental  SEA  results.  Before 
summarizing  these  results,  let  us  note  that  our  problem  formulation  thus  far  in  terms  of  a  Lyapunov 
equation  as  in  (1.6)  already  represents  a  point  of  departure  from  the  techniques  employed  in  [1-16]. 
Motivated  by  the  literature  on  large  scale  systems  theory  [39],  we  utilize  Kronecker  matrix  algebra 
[40,41]  and  M-matrix  theory  [39,42]  as  our  principal  mathematical  tools.  In  an  earlier  paper  [43]  we 
used  similar  tools  to  analyze  the  stability  and  performance  robustness  of  interconnected  systems. 
The  results  of  [43],  which  were  themselves  motivated  by  SEA,  thus  served  as  a  precursor  to  the 
SEA  extensions  given  in  the  present  paper. 

Perhaps  the  most  fundamental  tenet  of  SEA  is  that  quantities  such  as  IE[|y|3]  cab  be  estimated 
or  approximately  determined  solely  in  terms  of  the  “modal  energies”.  In  our  notation  the  modal 
energies  translate  into  the  real,  nonnegative  diagonal  elements 

£fc±gfc*=  1-..  IE[|zfc(3]  (1.9) 

*—♦00 

of  the  second-moment  matrix  Q.  For  example,  if  the  system  is  a  set  of  mechanical  subsystems  me¬ 
chanically  coupled,  then  E *.  corresponds  to  the  kinetic  or  potential  energy  of  one  of  the  vibrational 
modes  of  a  subsystem.  In  Section  2  we  discuss  the  various  conditions  under  which  it  suffices  to 
determine  the  Ek  in  order  to  evaluate  the  mean-square  response  of  quantities  of  interest. 

Having  argued  that  mean-square  response  measures  of  interest  can  be  deduced  from  knowledge 
of  the  modal  energies,  a  second  central  tenet  of  SEA  is  that  it  is  possible  to  formulate  a  set  of 
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n  linear  equations  that  involve  only  the  quantities  Ek  and  that  are  sufficient  to  determine  these 
quantities.  Note  that  the  diagonal  elements  of  (1.6)  give  power  flow  relations  of  the  form 


(2i>k  -  2Re  Hk)Ek 

+  Ilk 

Vkk, 

power  dissipated 

power  flow 

power  input 

by  the  Arth  mode 

from  the  Jbth  mode 

due  to  external 

due  to  damping 

to  all  other  modes 
due  to  coupling 

disturbances 

where  Ilk  is  given  by 

n 

Ilk  —  -  ^( GktQtk  +  QktGkt)- 

(-1 

l+h 

Statistical  Energy  Analysis  asserts  that  the  quantities  77*  can  be  evaluated  as  linear  functions  of 
the  Ei  s  alone  so  that  a  relationship  holds  of  the  form 

n 

iik  =  Pm € 

tel 


Thus,  if  (1.11)  holds  then  by  using  (1.10)  one  need  only  solve  n  linear  equations  for  the  modal 
energies  in  place  of  solving  the  ^n(n  +  1)  equations  corresponding  to  the  n  x  n  Lyapunov  equa¬ 
tion  (1.6).  Relation  (1.11)  has  been  demonstrated  in  several  special  cases,  namely,  two  coupled 
oscillators,  n  identical  oscillators  with  identical  coupling,  and  n  nonidentical  oscillators  with  weak 
inter-modal  coupling  [1-6].  In  Section  3,  however,  and  without  restrictions  (1.7),  (1.8),  we  use 
Kronecker  algebra  to  deduce  directly  from  (1.6)  that  the  modal  energies  Ek  are  determined  by  an 
energy  equation  of  the  form 

t n+P)E  =  V ,  (1.12) 


where 

H  =  diag  {nk),  Hk  —  2vk  -  2Re  Hk,  P  €  IRnxn, 


-vr 

rv„i 

r^ii 

Q  ii 

A. 

. 

<lll 

. 

. 

\ 

I 

* 

— 

* 

.Vn. 

-Vnn. 

.En. 

_Qnn  _ 

By  comparir*  (1.10)  to  (1.12)  it  can  be  seen  that  the  expression  (1.11)  for  Ilk  is  precisely  the  fcth 
element  of  PE.  So  long  as  the  overall  system  is  asymptotically  stable,  relations  of  the  form  (1.10) 
and  (1.11)  hold  regardless  of  the  number  of  modes  or  the  magnitude  of  the  couplings. 


Further  conditions  on  n  +  P  that  arise  in  the  cases  of  two  oscillators,  n  identical  oscillators 
with  identical  coupling,  or  nonidentical  oscillators  with  weak  coupling  lead  to  an  energy  difference 
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power  flow  proportionality  as  in  [11,12].  Specifically,  suppose  that 


Pkt<  0,  k,i  =  1, . . .  ,n,  (1.13) 


and 

Pkk  =  ^i\Pkt\,  1»  •••.*»•  (1-14) 

tel 

Then,  defining  <Jkt  =  |  ?kt\,  k  ^  t,  so  that  <Jkt  >  0,  it  follows  from  (1.11)  that 


nk  =  ^2<TktiEk-  Et).  (1-15) 

im  1 

L+h 

In  other  words,  power  flow  from  the  Jfcth  mode  to  all  other  modes  is  the  sum  of  the  individual 
power  flows  from  mode  k  to  mode  l,  which  are  proportional  to  the  energy  differences  Ek  -  Et- 
Note  that  power  always  flows  from  more  energetic  modes  to  less  energetic  modes  (because  of  the 
nonnegativity  of  the  coefficients  <Tkt)-  Substituting  (1.15)  into  (1.12)  yields 


-  £*)  =  V*,  (1.16) 

tml 

t+h 

which  is  an  energy  balance  relation.  Equations  (1.10)  and  (1.16),  which  govern  energy  exchange 
among  coupled  oscillators,  are  completely  analogous  to  the  equations  of  thermal  transfer  with  the 
modal  energies  playing  the  role  of  temperatures. 


In  physical  situations  involving  nonconservative  couplings,  we  show  that  although  (1.14)  no 
longer  holds,  it  is  still  possible  in  the  case  of  weak  couplings  to  obtain  a  generalized  power  flow 
proportionality.  In  this  case  there  exists  a  set  of  positive  scale  factors  Dk  >0,  k  =  1, . . . ,  n,  such 
that,  with  Ek  —  "p^Ek,  the  energy  difference  power  flow  proportionality  is  given  by 


nh  =  Y,*ki(Ek-Et),  (1.17) 

£-1 

t+h 


where  akt  —  Dt<*kt-  Note  that  (1.17)  is  not  merely  a  rewriting  of  (1.15)  since  in  general  Dk  #  Dt. 
With  (1.17),  the  energy  equation  (1.12)  assumes  the  form  of  a  generalized  energy  balance  relation 
given  by 

n 

A  kEk  +  ^  d'ktjEk  —  Et)  =  Vfc,  (1.18) 

te* 

where  k  =  l,...,n.  That  is,  there  is  a  set  of  re-scaled  energies  such  that  (1.12)  looks  like  the 
equations  of  thermal  transfer.  This  result,  given  in  Section  4,  generalizes  (1.15),  (1.16)  to  the  case 
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of  weak  but  otherwise  arbitrary  (not  necessarily  conservative)  modal  couplings.  These  results  are 
obtained  by  means  of  M-matrix  theory  [39,42]. 

While  deriving  energy  difference  power  flow  proportionality  relations,  we  show  that  the  explicit 
expressions  given  for  P  in  the  SEA  literature  are  actually  first-term  approximations  in  a  series 
expansion  for  P.  Indeed,  it  turns  out  that  P,  which  is  given  by  a  complicated  expression  involving 
t /,n,H,  and  G,  agrees  with  the  customary  SEA  expressions  for  “small”  G.  This  in  done  by  obtaining 
explicit  expressions  for  the  terms  of  a  series  expansion  of  P  in  ascending  powers  of  the  matrix 
elements  of  G. 

Since  the  modal  energies  satisfy  equations  analogous  to  those  of  thermal  transfer,  it  might  be 
expected  that  if  the  coupling  coefficients  Gkt  are  large  compared  to  the  modal  dampings,  then  the 
energies  should  be  approximately  equal,  that  is, 

Ei  as  Ej  —  •  •  ■  —  En.  (1.19) 

Section  7  provides  a  formulation  and  proof  of  this  “energy  equipartitioning”  phenomenon. 

At  this  point,  it  is  evident  that  this  paper  deals  only  with  certain  deterministic  aspects  of  SEA. 
Rigorous  exploration  and  extension  of  the  “Statistical”  aspect  of  Statistical  Energy  Analysis,  which 
addresses  the  possibility  of  uncertainties  in  the  system  parameters  and  coupling  coefficients,  will 
form  the  subject  of  a  future  paper.  Other  extensions  of  the  present  paper  are  briefly  mentioned  in 
Section  8. 
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2.  Characterization  of  System  Response  in  Terms  of  the  Modal  Energies 

Here  we  examine  the  conditions  under  which  it  suffices  to  compute  only  the  modal  energies 
(1.9)  in  order  to  estimate  response  quantities  such  as  limt_<»  IE[|y|3].  To  carry  out  the  necessary 
calculations,  we  shall  utilize  a  somewhat  unconventional  notation  for  the  diagonal  and  off-diagonal 
portions  of  a  matrix.  Specifically,  for  M  €  Cnxn  define 

{M}=  diag  (Affcfc),  (M)  =  M-  {M}. 

k= 1 . n 

For  convenience,  several  identities  involving  these  definitions  are  given  in  Appendix  A. 

Next  we  define  the  matrix 

A  =  -v+jn  +  H 

and  note  that 

A=  diag  (A*),  (2.1) 

fc=l . n 

where 

Ak  =  -Vk  +ffik  +  Hk- 
Then  the  Lyapunov  equation  (1.6)  becomes 

0  =  AQ  +  QA*  +GQ  +  QG*  +V.  (2.2) 

Using  the  identities  of  Appendix  A  to  decompose  the  Lyapunov  equation  (1.6)  into  its  diagonal 
and  off-diagonal  parts,  we  obtain,  (noting  A  =  (A)  and  G  =  (G)) 

0  =  A{Q}  +  {Q}A*  +{G{Q)  +  (Q)G*>  +  {V},  (2.3) 

0  =  A(Q)  +  {Q)A*  +  (G(Q)  +  (Q)G*>  +  G{Q>  +  (Q)G*  +  £0,  (2.4) 

while  (1.5)  becomes 

Um  lEdvl1!  =  tr[{C-CHQ>l  +  tr[(C-Q(Q>|.  (2.5) 

The  underlined  terms  in  (2.4)  and  (2.5)  are  zero  when  V  and  G*C  are  diagonal.  Furthermore,  they 
can  be  neglected  when  the  following  conditions  hold  either  separately  or  in  combination: 

t)  The  term  (V)  in  (2.4)  can  be  neglected  when  the  modal  excitation  forces  are  uncorrelated,  in 
which  case  (V)  ~  0.  This  occurs  when  excitations  are  spatially  distributed  with  very  short 
correlation  length. 
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it)  The  underlined  terms  in  (2.4)  and  (2.5)  are  negligible  when  Vk,vt  <<|  G *  -  Gt  |,  that  is,  the 
case  of  large  modal  frequency  separation  relative  to  modal  damping. 

Hi)  The  underlined  terms  in  (2.4)  and  (2.5)  are  negligible  in  the  case  of  a  distributed  structural 

system  with  very  high  modal  density  wherein,  for  fixed  k  and  for  l  =  1 . n,  the  real  and 

imaginary  parts  of  Gu,Gtk,  and  Vki  have  many  sign  reversals  for  modes  within  any  narrow 
frequency  band.  These  sign  reversals  essentially  cancel  out  the  contributions  of  (Q)  in  (2.5) 
and  the  effect  of  (V)  on  {<?}  in  (2.4)  (see  [13]  for  details). 

When  conditions  *)-ro)  are  satisfied  either  separately  or  in  combination,  we  have  the  approxi¬ 
mate  equations 


0  =  A{Q }  +  {Q}A*  +  {G(Q)  +  (< Q)G *}  +  {V},  (2.6) 

0  =  A{Q)  +  <Q)A*  +  (G(Q)  +  (Q)Gm)  +  G{Q>  +  {Q}<T ,  (2.7) 

lim  IE[|y|J]  =  tr[{C*C}{<3}].  (2.8) 

t— *oo 

These  equations  are  exact  when  V  and  C*C  are  diagonal;  they  are  good  approximations  under 
conditions  i)-tii).  Note  that  these  approximations  have  no  impact  on  stability  analysis. 

The  salient  feature  of  (2.8)  is  that  the  response  quantity  involving  y  can  be  expressed  in  terms 
of  the  diagonal  elements  of  Q.  As  mentioned  in  Section  1,  the  diagonal  elements  Qkk  have  the 
physical  significance  of  either  kinetic  or  potential  energies  of  the  vibrational  modes.  Although  we 
need  only  calculate  the  n  diagonal  elements  of  Q  (the  “system  modal  energies”)  to  evaluate  (2.8), 
it  is  still  apparently  necessary  to  solve  an  n  x  n  Lyapunov  equation  to  obtain  all  of  Q.  In  fact, 
however,  we  now  proceed  to  use  Kronecker  matrix  algebra  to  eliminate  the  off-diagonal  part  (Q) 
from  (2.6)  and  (2.7),  thereby  producing  a  system  of  only  n  equations  determining  {<?},  rather  than 
the  \n(n  +  1)  equations  that  characterize  all  of  Q. 
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3.  Determination  of  the  Modal  Energy  Equations 


Here  we  show  that  the  decomposed  Lyapunov  equation  (2.6),  (2.7)  can  be  reduced  to  a  system 
of  n  equations  involving  only  the  modal  energies  Ek,  that  is,  the  diagonal  elements  of  Q.  To  do 
this  we  employ  the  Kronecker  matrix  algebra,  the  basic  definitions  and  identities  of  which  are 
summarized  in  Appendix  B.  Note  that  the  basic  operators  are  the  vee  operator,  which  stacks  the 
columns  of  a  matrix  into  a  vector,  and  the  vecd  operator,  which  stacks  only  the  diagonal  entries  of  a 
square  matrix  into  a  vector.  Appendix  B  reviews  the  definition  of  the  Kronecker  product  and  sum 
along  with  identities  (B.l)  through  (B.10),  which  are  well  known  [40,41].  The  remaining  identities 
(B.11)-(B.18)  are  new  and  their  proof  is  left  to  the  reader.  The  matrices  £  and  £±  are  diagonal 
projection  matrices  that  allow  us  to  separate  the  entries  of  vec  M  corresponding  to  the  diagonal 
and  off-diagonal  elements  of  a  square  matrix  M  (see  (B.13)  and  (B.14)).  We  can  now  define 


r^ii 

’  Qn ' 

En 

=  vecd  Q  — 

Qt  a 

.En. 

rv.i  i 

•  Qnn  - 

■Vul 

V  £  vecd  V  = 


Vt 

i 


VM 

Vnn 


(3.1) 

(3.2) 


which  are  real,  nonnegative  vectors  since  Q  and  V  are  Hermitian  nonnegative-definite  matrices. 
Furthermore,  define 


diag  (nk),  /ifc  =  2i/fc  -  2Re  Hfc.  (3.3) 

Jbsl,...,n 

and  note  that  /j  =  —(A  +  A)  =  — 2Re  A. 


Theorem  3.1.  Assume  that  A  +  G  is  asymptotically  stable,  let  Q  €  <DnXn  be  the  unique 
Hermitian  nonnegative-definite  solution  to  the  Lyapunov  equation  (1.6),  and  define  the  nonnegative 
vectors  E,  V  €  IR"  by  (3.1),  (3.2).  If  A  ©  A  +  £±(6  ©  G)£±  is  nonsingular,  then  E  and  V  satisfy 

0 i  +  P)E  =  V,  (3.4) 

where  m  €  IRnxn  is  defined  by  (3.3)  and  P  6  IRnxn  is  defined  by 

P  =  iT{G  ®G)£±[A®  A+  £±{G  ®G)£X}~1£±(G  ®G)i .  (3.5) 

Furthermore,  n  +  P  is  nonsingular  and  its  inverse  (fi  +  P)~1  is  a  nonnegative  matrix.  Finally, 

n 

lim  IE[|y|a]  =  £  \Ck\2Ek.  (3.6) 

**♦00  1 

k=  1 
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(3.7) 


Proof.  Applying  the  vec  operator  to  (2.6)  and  using  (B.7),  (B.13),  and  (B.16)  yields 
0  =  (A  ©  A)vec{Q}  +  vec {G{Q)  +  (Q)G*}  4-  vec{^  } 

=  (A  ©  A)/ vecd  Q  +  £  (<5  ©  G)vec(g)  +  £vecd  V. 

Next  applying  the  vec  operator  to  (2.7)  and  using  (B.7),  (B.14),  (B.15),  (B.16),  and  £x(A©  A)  = 
(A  ©  A)<fx,  yields 

0  =  (A  ©  A)vec(Q)  4-  vec(G(<3)  4-  (Q)G*)  +  (G  ©  G)vec{g} 

=  (A  ©  A)vec(Q)  4-  £x(G  ©  G)£±vec(g)  4-  (G  ©  G)l vecd  Q  (3.8) 

=  [A©  A  4-  £X{G  ©  G)£x]vec(Q)  4-  £X(G  ©  G)£vecd  Q. 

Since  A  ©  A  4-  £±.(G  ©  G)<?x  is  assumed  to  be  nonsingular,  (3.8)  and  (B.14)  imply 

vec (Q)  =  — £x(A©  A  4-  £X(G  ©  G)£xj_1£x(6  ©  G)£vecd  Q.  (3.9) 

Substituting  (3.9)  into  (3.7)  yields 

0  =  [(A  ©  A)i  -£{G®  G)£X[A  ©  A  +  £±(G  ©  G)£x\~ltx{G  ©  G)£ jvecd  Q  + l  vecd  V.  (3.10) 
Next  note  that 

£t(A  ©  A)i  =  A  4-  A  =  — /i.  (3.11) 

Multiplying  (3.10)  by  £ T  and  using  (3.11),  (B.17),  and  (B.18)  yields 

0  =  -[m  4-  iT(G  ©  G)£x[A  ©  A  +  £±{G  ©  G)£X]~1£X{G  ©  G)£ jvecd  Q  4-  vecd  V, 


or,  using  (3.5), 


which  is  (3.4). 


(n  4-  P  Jvecd  Q  —  vecd  V, 


To  show  that  the  n  x  n  matrix  P  defined  by  (3.5)  is  real,  take  the  complex  conjugate  of  (3.5) 
and  use  (B.8)  and  (B.10)-(B.12)  to  obtain 

P  =  £t(G  ©  G)£X[A  ©  A  4*  £x[G  ©  @)£x]~l£x(G  ©  G)£ 

=  iTU{G  ©  G)U£±[U{A  ©  A)U  +  £x{0  ©  G)U£X]~1£XU{G  ©  G)U£ 

—  £t(G  ©  G)fx[A  ©  A  4"  £x{G  ©  G)£x\  1  Ex{G  ©  G)£ 

=  P. 
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Next,  to  show  that  ft  +  P  is  nonsingular,  note  that  since  V  €  IRn  is  an  arbitrary  nonnegative 
vector,  the  rank  of  n  +  P  is  n.  Thus  n  +  P  is  nonsingular.  Furthermore,  it  can  be  seen  that  if  V 
is  the  sth  column  of  /n>  then  the  nonnegative  solution  E  of  (3.4)  is  the  tth  column  of  (p  +  P)~l. 
Hence  (/t  +  P)~l  is  a  nonnegative  matrix.  Finally,  (3.6)  follows  from  (2.8).  Q 

Remark  3.1.  Suppose  that  G  is  symmetric,  that  is,  G  =  GT ,  but  not  necessarily  real.  Then 
using  (B.3)  it  is  easy  to  show  that  P  (which  is  real)  is  also  symmetric.  Hence  in  this  case  n  +  P 
and  (/i  +  P)~l  are  both  real  symmetric  matrices. 


As  a  separate  result  we  state  the  following  converse  of  Theorem  3.1. 

Proposition  3.1.  Assume  that  A®  A+£_l(<5®G)£_l  is  nonsingular,  let  V,n,  and  P  be  defined 


by  (3.2),  (3.3),  and  (3.5),  and  suppose  there  exists  a  nonnegative  solution  E  = 


e  nt*  to 


lEnl 


equation  (3.4).  Then  the  matrix  Q  €  €nxn  defined  by 


Q=  diag  (fi’<)  +  vec-1((A©A  +  ^x(<5©G)£±]-1fx(G©G)££:),  (3.12) 

»=i . » 

is  Hermitian  and  satisfies  (1.6).  If,  in  addition,  Vjb  >  0,  k  =  1 . .  and  Q  is  positive  definite, 

then  A  +  G  is  asymptotically  stable. 

Proof.  The  fact  that  Q  given  by  (3.12)  satisfies  (1.6)  follows  by  reversing  the  algebraic  steps 
leading  to  (3.4).  To  show  that  Q  is  Hermitian,  note  that  using  (B.8)-(B.10)  we  have 

Q*=  diag  (£?<)+ [vec-l(fx(A©A+ fx(<3©G)fx]_1fx (G©5)£* 5)] T 

„n 

=  diag  (Ei)+[vec-l(£±[U{A®A)U  +  £±U{G®G)U£±}-l£xU(G®G)UiE)]T 

»=X,...,n 

=  diag  (£<)+  [vec-1(£x^[A©A  +  £x(<5©G)fx]"1^x(<5©G)£F;)]T 

. *» 

=  diag  (Ei)+  \ec~l (£x[A  ®  A+  £±(G  ®G)£1]~1£±(€r  QG)i E) 

. n 

=  Q 

Finally,  the  stability  of  A  +  G  follows  from  standard  Lyapunov  theory  [44,  Lemma  12.2].  □ 

Theorem  3.1  and  Proposition  3.1  show  that  for  the  purpose  of  determining  the  diagonal  entries 
of  Q,  that  is,  the  modal  energies  Ei,.. . ,  E„,  equation  (3.4)  is  equivalent  to  equations  (2.6)  and 
(2.7).  This  verifies  the  tenet  of  Statistical  Energy  Analysis  that  there  exists  a  system  of  n  linear 
equations  that  determine  the  modal  energies  alone.  Moreover,  comparing  the  fcth  equation  in  (3.4) 
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with  the  power  balance  relation  (1.10),  we  now  see  that  Ilk,  the  power  flow  from  the  fcth  mode  to 
all  other  modes  due  to  coupling,  is  given  by 

n 

nk  =  J2  PkiEt,  (3.13) 

£=1 

where  P  €  IR  is  the  {k,t)  element  of  P.  Thus  the  expression  (1.11)  is  also  verified.  In  the  next 
section  we  further  explore  the  structure  of  Ilk  and  P  to  derive  a  generalization  of  the  energy 
difference  power  flow  proportionality  (1.15)  for  weak  but  arbitrary  coupling  matrices  G. 
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4.  Analysis  of  the  Energy  Equation:  Energy  Difference  Power  Flow 

Proportionality 

In  this  section  we  analyze  the  energy  equation  (1-12)  to  determine  conditions  under  which  an 
energy  difference  power  Sow  proportionality  holds.  Under  the  assumption  that  the  off-diagonal 
elements  of  P  are  nonpositive,  we  obtain  a  generalized  power  Sow  proportionality  involving  scaled 
model  energies.  In  Section  5  we  then  show  that  this  result  holds  for  weak,  but  otherwise  arbitrary, 
couplings.  Specializing  further  in  Section  6  to  the  conservative  case  involving  skew-Hermitian 
couplings,  we  obtain  a  power  Sow  proportionality  involving  the  actual  (unsealed)  modal  energies. 

The  development  requires  several  deSnitions  and  results  from  matrix  theory  [39,42].  A  matrix 
M  €  IRnxn  is  called  a  Z-matrix  if  Mkt  <0,  k  ^  t  k,t  —  1, . . .  ,n.  Note  that  a  Z-matrix 
M  €  IRnxn  can  always  be  placed  in  the  form 


M  =  al-N,  (4.1) 

where  a  >  0  and  N  >>  0,  N  G  IRnxn.  If  (4.1)  can  be  satisSed  with  a  >  p{N)  (p  denotes  spectral 
radius),  then  Af  is  called  an  M -matrix.  If,  furthermore,  a  >  p(N),  then,  since  det  M  0,  M  is 
a  nonsingular  M-matrix.  There  are  numerous  (at  least  50)  equivalent  conditions  under  which  a 
Z-matrix  is  a  nonsingular  M-matrix  [42].  We  now  summarize  those  conditions  that  will  be  used 
here.  We  shall  call  B  G  IRnxn  diagonally  dominant  if 

n 

fc  =  l,...,n.  (4.2) 

\ 

t+h 


Lemma  4.1.  Let  M  G  IRnxn  be  a  Z-matrix.  Then  the  following  are  equivalent: 
t)  Af  is  a  nonsingular  M-matrix, 
it)  M  is  nonsingular  and  M-1  >>  0, 
tit)  the  real  part  of  each  eigenvalue  of  Af  is  positive, 

tv)  there  exists  positive  diagonal  P  €  IRnxn  such  that  MD  is  diagonally  dominant. 

Proof.  See  conditions  {N39),  (G20),  and  (M35)  on  pages  134-138  of  [42].  □ 

Returning  to  the  energy  equation  (1.12),  we  focus  on  the  coefficient  matrix  p  +  P .  The  crucial 
condition  that  p+  ?  is  a  Z-matrix  will  be  shown  later  for  the  case  of  weak,  but  otherwise  arbitrary, 
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couplings.  First  we  recall  from  Theorem  3.1  that,  under  the  assumptions  of  that  Theorem,  n  +  P 
is  nonsingular  and  (p  +  P)~1  >>  0.  Thus,  condition  ii)  of  Lemma  4.1  with  M  =  n+  P  can  be 
invoked  to  yield  conditions  t),  tit)  and  tv). 

Proposition  4.1.  Suppose  that  the  assumptions  of  Theorem  3.1  are  satisfied  and  assume  that 
P  is  a  Z-matrix.  Then 

t)  n  +  P  is  a  nonsingular  M-matrix, 

ii)  the  reed  part  of  each  eigenvalue  of  p  4-  P  is  positive, 

tv)  there  exists  positive  scalars  Di , . . . ,  Dn  such  that 

n 

Dk(fik  + Pkk)  >^>2Dl  |^m|,  fc  =  1, ...,n.  (4.3) 

<■1 

Proof.  First  note  that  since  n  is  a  diagonal  matrix,  n  +  P  is  a  Z-matrix  if  and  only  if  P  is 
a  Z-matrix.  Since,  by  Theorem  3.1,  (/*  +  P)~l  >>  0,  condition  it)  of  Lemma  4.1  is  satisfied  with 
M  =  fi  +  P.  Hence  conditions  t),  Hi),  and  iv)  of  Lemma  4.1  are  also  satisfied.  Now  it  need  only  be 
noted  that  (4.3)  is  equivalent  to  (4.2)  with  B  =  (*t+  P)D  and  D  —  diagfc_x . n(Dfc).  □ 

Remark  4.1.  Suppose  G  is  symmetric  but  not  necessarily  real.  Then  by  Remark  3.1  P  is 
symmetric.  Since  ix  +  P  is  also  symmetric,  fi  +  P  has  only  real  eigenvalues.  It  thus  follows  from 
condition  it)  of  Proposition  4.1  that  n  +  P  has  only  real  positive  eigenvalues.  Hence  in  this  case 
H  +  P  is  a  symmetric  positive-definite  matrix. 

Corollary  4.1.  Suppose  that  the  assumptions  of  Theorem  3.1  are  satisfied  and  assume  that 
P  is  a  Z-matrix.  Then  there  exist  positive  scalars  fik  >  0,  k  =  1, . . .  ,n,  and  nonnegative  scalars 
°ki  >0,  k  ^  t ,  k,  l  =  1, . . . ,  n,  such  that 

n 

fikEk  +  ^ &ki{Ek  -  Et)  =  V b,  k=l,...,n,  (4.4) 

tm  1 

where  Ek  =  -£;Ek,  k  =  l,...,n. 

Proof.  Using  (4.3)  of  Proposition  4.1,  define  flk  >  0  by 

n 

Afc  —  Dk(fik  4-  Pkk)  -  ^  Dt  | Pki\y  k  =  l,...,n.  (4.5) 

t+k 
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Next  note  that  with  Ek  =  -^Ek  and,  since  P  is  assumed  to  be  a  Z-matrix,  Pkt  =  -  \  Pkt\,  k  ^  i, 
the  &th  equation  of  (1-12)  yields 

n 

Dk[f*h  +  Pkk)Ek  ~  TM  Pki\ Et  =  Vk-  (4.6) 

Iml 

l+k 


Combining  (4.5)  and  (4.6)  yields 


fikEk  +  Y^Dt\Pkt\(Ek  -  Et)  =  V*,  (4.7) 

<■1 


which  implies  (4.4)  with  arkt  =  Dt  \  Pkt  |  •  □ 


Equation  (4.4)  can  be  viewed  as  a  generalized  energy  balance  relation  since  it  involves  scaled 
modal  energies  rather  than  the  modal  energies  themselves.  Furthermore,  comparing  (4.4)  to  (1.10) 
it  follows  that 

n 

Ilk  =  (pk  ~  Dknk)Ek  +  Oki{Ek  -  Et).  (4.8) 

tml 

l+h 

That  is,  the  power  flow  from  the  fcth  mode  to  all  other  modes  is,  aside  from  the  offset  term 
{tik  -  Dknk)Ek ,  proportional  to  the  difference  between  scaled  modal  energies.  In  Section  6  we  show 
that  under  conservative  couplings  (4.8)  becomes  an  actual  (nonscaled)  energy  difference  power  flow 
proportionality.  Next,  however,  we  show  that  P  is  a  Z-matrix  under  weak  coupling. 
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5.  Analysis  of  P  in  the  Case  of  Weak  Coupling 

In  the  previous  section  the  generalized  power  flow  proportionality  was  based  upon  the  assump¬ 
tion  that  P  is  a  Z-matrix.  In  this  section  we  show  that  this  assumption  is  valid  in  the  case  of  weak 
but  otherwise  arbitrary  couplings.  To  do  this  we  expand  P  in  terms  of  powers  of  G  and  then  show 
that  the  first  term  in  the  expansion  is  a  Z-matrix. 

To  begin  we  define  for  convenience 

A  =  A®A,  9  =  G®G  (5.1) 

so  that  P  defined  by  (3.5)  can  be  written  as 

P  =  iT9£x(A  +  £x9S±)-l£x9i-  (5-2) 

For  r  =  0, 1, 2, . . . ,  it  is  easy  to  confirm  the  identity 

(A  +  £x9£±r1  =  Y,{-A-l£x9Z±)i*-x  +  (-A-1£x9£x)r+1(A  +  £x9^)~l.  (5.3) 

»=0 

Combining  (5.2)  and  (5.3)  it  follows  that 


i=0 

(5.4) 

where 

Pi  =  iT9£x(-A-l£x9£±Y*~1£x9i 

(5.5) 

and 

Zr  =  ir9£x(-A~1£x9£. L)r+1(X  +  £x9£x)~l£x9t- 

(5.6) 

Note  that 

n«rii=o(iM-ifxprxirs) 

(5.7) 

for  ||p ||  — »  0  (where  ||  •  |j  denotes  arbitrary  submultiplicative  matrix  norms), 
incurred  in  approximating  P  by  Pi  depends  on  the  size  of  A~1£±$£x- 

with  the  SEA  literature  since  A~1£x9ix  can  be  viewed  as  the  ratio  of  modal 

damping.  For  *=0,1  we  have 

Clearly  the  error 

This  is  consistent 

coupling  to  modal 

Po  =  iT9£xA~1£x5£, 

Pi  =  -( !’T9£xA-l£x9£xA~l£x9£- 

(5.8) 

(5.9) 
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Lemma  5.1.  Suppose  that  the  assumptions  of  Theorem  3.1  are  satisfied.  Then 

n 

(Po)kk  =  -2Re[y^Cfct<5t*A*]i  fc=l,...,n,  (5.10) 

t=i 

(Po)ki  =  -2\Gu\*Re  TtM,  l?k,  £, *  =  1, . . .  ,n,  (5.11) 

n 

=  -2Re[  ^  C?wAjfcG!^mrm*Gmjfc],  fc=l,...,n,  (5-12) 

£,m= 1 

n 

(Pl)kt  =  -2Re  [  ^  (G km.rmkG mirtkG kl  +  G  km,Fm.kFmiGmtG  kt> 

m= 1 

-H  ^fciAfc  AmGmiGfcm)]  ,  £  7^  fc,  £,  k  =  l,...,n,  (5.13) 

where 

A*  =  ■  j  =  [*'£  +  *'*  -  {Ht  +  Hk)  +}(Gk  —  Gt)]~l.  (514) 

Proposition  5.1.  Suppose  that  the  assumptions  of  Theorem  3.1  are  satisfied  and  assume, 
furthermore,  that  A**  +  fit  >  0,  k  ^  l,  k,l  =  1, . . . ,  n.  Then  P0  is  a  Z-matrix.  If,  furthermore, 
+  k  l,  fc,/=l,...,n,  Gw  ^  0,  fc  =£  £,  k,l  =  1, ...  ,n,  and  ||G||  is  sufficiently 

small,  then  P  is  a  Z-matrix. 

Proof.  From  (5.11)  we  have 

(Po)kt  —  -2 1  Gkt\3  Re  A* 

—  ~  I  Gkt  |J  0*fc  +  Mi)/  I  A*  |S 

<  0. 

Hence,  P0  is  a  Z-matrix.  If,  in  addition,  Mfc  +  Pt  >  0  and  Gkt  ^  0,  then  ( Po)ki  <  0.  In  this  case 
||G||  sufficiently  small  implies  Pkt  <  0  so  that  P  is  a  Z-matrix.  □ 

To  understand  the  significance  of  Proposition  5.1,  consider  in  place  of  (3.4)  the  approximate 
energy  equation 

(m  +  P0)E  =  V.  (5.15) 

If  ||G||  is  small,  that  is,  the  coupling  G  is  weak,  then  the  norm  of  the  residual  Pq  is  of  order 
IM-1£x5 A.II3-  Hence  in  this  case  (5.15)  can  serve  as  an  approximation  to  (3.4). 

Proposition  5.1  also  shows  that  if  fik  +  Hi  >  0  and  Gkt  #  0 ,  k  ^  l,  k,  l  =  1, . . . ,  n,  then  P 
itself  is  a  Z-matrix  so  that  Corollary  4.1  can  be  applied.  Note  that  if  Gkt  —  0  then  (Po)kt  =  0  and 
thus  the  sign  of  Pkt  depends  on  higher  order  terms  in  the  expansion  of  Pkt ■  It  is  interesting  to  note 
that  if  Gkt  —  0  then  (Pi)kt  is  also  zero  so  that  in  this  case  terms  of  even  higher  order  play  a  role. 
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6.  Specialization  to  the  Case  of  Conservative  Couplings 

Many  physical  situations  involve  only  passive  or  energy-conservative  couplings  among  subsys¬ 
tems.  This  is  the  case  considered  in  the  SEA  literature.  To  model  this  situation  we  assume  that  H 
and  G  are  skew-Hermitian.  If  V(x)  =  x*x  represents  the  total  energy  of  the  system,  then  it  follows 
that  energy  dissipation  along  trajectories  of  the  system  (1.3)  with  w  =  0  is  given  by 

=  -2x*vx  <0,  x  ^  0,  (6.1) 

at 

which  is  identical  to  the  energy  dissipation  of  the  uncoupled  system.  Thus  skew-Hermitian  coupling 
has  no  effect  on  the  total  system  energy.  To  analyze  this  case  we  begin  with  the  following  lemma 
which  corresponds  to  equation  (6)  of  [ll]. 

Lemma  6.1.  Suppose  that  the  assumptions  of  Theorem  3.1  are  satisfied  and,  furthermore, 


assume  that  G  is  skew-Hermitian.  Then 


Pe  =  0, 


where  e  =  .  I  . 


Proof.  It  suffices  to  show  that  (G  ©  G)£e  =  0.  Note 

(G  ©  G)le  =  (G  ©  G)i vecd  In 
=  (G  ©  G)vec{/„} 

=  (G  ©  G)vec  In 
=  vec(G*  +G) 

=  0.  □ 

Since  e  ^  0,  P  has  a  nontrivial  nullspace  and  thus  (6.2)  implies  that  P  is  singular.  Note  that 
(6.2)  can  be  written  as 

n 

?kt  =  °,  fc  =  1, . . . ,  n.  (6.3) 

4=1 

If,  in  addition,  P  is  a  Z-matrix,  then  (6.3)  is  equivalent  to 


Pkk  =  Y,\Pki\,  k=  l,...,n. 
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Defining  <Tkt  =|  Pkt\—  —  Pki,  k  ^  t,  k,l  =  1, . . .  ,n,  it  thus  follows  from  (1.11) 


Ilk  =  u(Ek  ~  Et)-  (6.5) 

<•1 

t+h 

Using  (6.5)  we  can  now  obtain  an  energy  difference  power  flow  proportionality  as  a  specialization 
of  (4.4).  This  results  is  obtained  directly  and  not  by  means  of  M-matrix  theory  which  was  used  to 
derive  (4.4). 

Proposition  6.1.  Suppose  that  the  assumptions  of  Theorem  3.1  are  satisfied,  assume  that  P 
is  a  Z-matrix,  and  that  G  is  skew-Hermitian.  Then  with  ffkt  =  | Pkt\ ,  k  ±  4,  k,l  =  1 . . .  ,n,  it 

follows  that 

n 

HkEk  +  ^^<rkt{Ek  -  Et)  —  Vk,  fc=l,...,n.  (6.6) 

4-1 

Proof.  Equation  (6.6)  is  the  fcth  equation  of  (3.4)  using  (6.4).  □ 

Remark  6.1.  Proposition  6.1  does  not  state  the  Hk  >  0,  which  is  needed  for  (6.6)  to  have  the 
physical  interpretation  of  an  energy  balance  relation.  Note,  however,  that  in  (4.4)  the  coefficient 
pLk  was  shown  to  be  positive  by  means  of  the  diagonal  dominance  characterization  of  nonsingular 
M-matrices.  Invoking  this  condition  here  would  lead  to  a  scaled  energy  balance  relation  in  place  of 
(6.7). 

Remark  6.2.  Suppose  in  addition  to  the  assumption  that  G  is  Skew-Hermitian,  we  assume 
that  Re  G  =  0.  Then  G  =  jG,  where  G  is  a  real  symmetric  matrix.  Consequently,  G  is  symmetric 
and  thus  Remark  3.1  implies  that  P  is  symmetric.  Hence  <7kt  =  otk  which  shows  that  the  power 
flow  from  the  kih  mode  to  the  4th  mode  is  equal  to  minus  the  power  flow  from  the  4th  mode  to  the 
fcth  mode. 
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7.  Equipartition  of  Energy 

In  the  case  of  conservative  couplings  as  considered  in  Section  6  we  can  show  that  energy 
equipartition  occurs  in  the  limit  of  strong  coupling.  Hence  assume  that  G  is  a  fixed  skew-Hermitian 
coupling  matrix  and  scale  G  by  7  >  0  so  that  (1.12)  is  replaced  by 

(M  +  P(7))£  =  V\  (7.1) 

where 

ph) = (7.2) 

We  are  interested  in  evaluating 

E=  lim  (11+ Ph))-1#.  (7.3) 

Tf-*oo 

We  sketch  the  main  steps  of  the  derivation. 

It  can  be  shown  that 

£=  lim  (li  +  iP)-1#.  (7.4) 

7-*  OO 

Now  assume  G  is  symmetric  so  that  P  is  symmetric.  Then  by  Corollary  7.6.3  of  [45]  it  follows  that 
bm  (p  +  7^)-1  =  p-1  -  p~l P n~$ P ,  (7.5) 

7— *oo 

where  (  )+  denotes  Moore-Penrose  generalized  inverse  of  P  (or  Drazin  generalized  inverse  since  P 

is  symmetric).  Now  suppose  that  G  is  also  skew- Hermitian .  Then  by  Lemma  6.1,  Pt  —  0  so  that 

P  pi  e  =  0  (7.6) 

Next  suppose  that  P  is  a  Z-matrix.  Then  it  follows  from  Lemma  6.4.1  of  [42]  (by  setting  fx  =  el) 
that  P  is  an  M-matrix.  Next  assume  P  is  irreducible,  which  is  the  case  if  all  modes  are  mutually 
coupled.  Then,  since  P  is  a  singular  irreducible  M-matrix,  it  follows  from  Theorem  6.4.16  of  [42] 
that  rank  P  =  n  —  1.  Thus  the  null  space  of  P  is  the  one-dimensional  subspace  spanned  by  e.  Now 
it  can  be  seen  that 


I  -  Pn~i(n~t  Pn~$)+  =  —  . 

(7.7) 

Hence  (7.4),  (7.5),  and  (7.6)  yield 

=  =  i, 

e/xe  efxe  'e^/ie7 

M 

Hence 

eTV 

—  £/2  —  ■  •  •  —  En  —  -  ■ 

eTiie 

(7.9) 

which  is  an  equipartition  of  energy. 
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8. 


Concluding  Remarks 


There  are  several  issues  and  questions  that  remain  to  be  explored: 

1.  How  restrictive  is  the  assumption  that  A©  A  +  £_l(G©G)£_l  is  nonsingular?  Can  the  inverse  of 
this  matrix  be  replaced  by  the  inverse  of  a  matrix  of  dimension  (n2  -  n)  x  (n2  —  n)  to  account 
for  the  rank  of  £±i 

2.  It  may  be  possible  to  redevelop  the  theory  with  real  (as  opposed  to  complex)  models  by  allowing 
nonscalar  blocks  in  the  interconnection  structure.  The  block  Kronecker  product  [46]  may  be 
useful  for  such  a  formulation. 

3.  Further  quantification  of  conditions  *)-***)  of  Section  2  may  be  useful.  The  theory  may  also  be 
extendable  to  the  case  (V)  ^  0. 

4.  It  may  be  possible  to  develop  transient  (as  opposed  to  steady-state)  results  for  power  flow. 

5.  It  is  well  known  that  power  flow  can  be  modeled  by  time-averaging  the  unforced  response  of 
the  system.  Such  a  dual  theory  may  provide  further  insights  into  the  power  flow  phenomenon. 
Note  that  a  time-averaging  theory  may  require  a  dynamic  model  that  is  conservative  rather 
than  asymptotically  stable. 

6.  Further  analysis  may  reveal  more  general  conditions  under  which  P  is  a  Z-matrix,  particularly 
for  the  case  of  strong  coupling. 

Acknowledgement.  We  wish  to  thank  Linda  Smith  for  transforming  the  original  manuscript 
of  this  paper  into  TfcX. 


22 


Appendix  A.  Identities  Involving  {•}  and  {•). 

For  matrices  A,  B  €  Cnx ,  the  following  identities  axe  satisfied: 

A  =  {A}  +  (A),  (A.l) 

A  =  {A}  O  (A)  =  0,  (A.2) 

A  =  (A)  O  {A}  =  0,  (A.3) 

{<A>}  =  0,  <{A»  =  0,  (A.4) 

{<A){i?}}  =  {{A}<B)}=0,  (A.5) 

{AB}  =  {{A}{B}  +  {A)(B}},  (A.6) 

((A){B}>  =  (A){B),  ({A)(B))  =  {A}(B),  {{A}{B})  =  0,  (A.7) 

<AB>  =  {A}(B)  +  (A)  {23}  +  ((A)(5)).  (A.8) 


Appendix  B.  Kronecker  Matrix  Algebra,  Definitions  and  Identities 
The  following  are  basic  definitions  and  identities: 
vec  and  vec-1  Operators:  For  A  €  Cnx"*, 

Aji 


An  i 
An 


vec  A  ==  I  An 


Ana 


llm 


,  vec“1(vec  A)  =  A 


vecd  Operator:  For  A  €  <Cnxn, 


vecd  A  = 


r  An 
Ajj 

l- 
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Kronecker  Product:  For  A  E  £nxm  and  B  E  <DpX<J, 


A  ®  B  = 


AuB  AijB  ...  AimB 
AaiB  AjjB  ...  AamB 

L  An  i  h  AnaB  ...  AnmB  J 


g  IRnpxm<l 


Kronecker  Sum:  For  A  E  CnXn  and  B  E  <DmXm, 


A©B  =  A®Jm  +  /n®B€  <nnmxnm 


Kronecker  Algebra  Identities:  For  compatible  complex  matrices  A,  B,C,D  : 


(A  +  B)®C  =  A®C  +  B&C,  (B.l) 

A®(B  +  C)  =  A&B  +  A&C,  (B.2) 

(A®B)T  =  AT®BT,  (A©B)t  =  At©Bt,  (B.  3) 

(A  ®  B)(C  ®  D)  =  (AC)  ®  (BD),  (BA) 

(A®  B)~l  =  A~l  ®  B~l,  (B.5) 

vec  ABC  =  (CT  ®  A)vec  B,  (S.6) 

vec (AB  +  BC)  =  (CT  ©  A)vec  B.  (B.7) 


Define  the  following  special  vectors  and  matrices  whose  dimensions  will  be  inferred  from  the  context 
in  which  they  are  used: 

er  =  column  vector  whose  rth  element  is  1  and  which  is  zero  otherwise, 

Er,  =  matrix  whose  (r,  a)-element  is  1  and  which  is  zero  otherwise  (Er,  —  erej), 

#  =  E,,E..88,'„ 

Er  =  Err,  where  Err  is  square, 

£=£r£r®e  r,  f=Er£r®£r,  i  JL  =  J  -  t . 

The  following  identities  hold  for  compatible  matrices  A,  B: 

U-1  =  ur  =  U,  (B.8) 
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vec  At  =  C/ve c  A, 

(B- 9) 

A  ®  £  =  U(B  ®  A)U,  A®  B  =  U[B  ©  A)U, 

(£.10) 

U£±U  =  £±t  U£±  =  £±U, 

(£.11) 

Pu  =  iTt  U£  =  £, 

(£.12) 

vec{A}  =  fvec  A  =  fvec {A}, 

(£.13) 

vec(A)  =  £±vec  A  =  £x,vec(A)i 

(£.14) 

£  =  £t  =  £2,  £j.  =  £l  =  £l 

(£.15) 

vec{A}  =  £vecd  A, 

(£.16) 

n = /, 

(£.17) 

ir£  =  r. 

(£.18) 
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1.  Background  and  Motivation 

Th«  nonlinear  compensator  deaign  introduced  in  Section  2  and 
subsequently  explored  in  the  remainder  of  the  paper  waa  initially 
motivated  by  the  problem  of  ayntheaiaing  control  algorithms  for 
vibration  suppression  in  large  flexible  structures.  Thus,  to  provide 
the  basic  background  for  the  present  development,  consider  a  flex¬ 
ible  structure  instrumented,  for  vibration  control  purposes,  with 
electromechanical  actuators  (to  provide  control  forces)  and  elec¬ 
tronic  sensors  (to  provide  measurements  of  structural  motion  used 
to  construct  appropriate  drive  signals  for  the  actuators).  In  terms 
of  modal  coordinates,  the  structural  plant  model  may  be  given  as: 

*-[i  *=(*0eEj"' ueKi(i) 

where 

n  =  diag  (17a)  =  modal  frequencies 
q  =  diag  (q*)  *=  modal  damping  ratios 

6  =  modal  actuator  influence  coefficients,  where  far  simplicity, 
we  consider  only  one  actuator  so  that  4  €  IR' 

Again,  for  simplicity  in  the  present  exposition,  we  suppose  that 
there  ia  one  rate  sensor,  collocated  with  the  actuator.  Then  the 
sensor  output,  y,  is  given  by: 

y  =  4Tx,  (2) 

The  actuator  input  signal,  u,  is  generally  synthesised  from  the 
measurement  signal,  y.  The  generic  form  of  a  linear  controller  is: 

u  =  -\Kp,  (3) 


The  physical  significance  of  £a  is  that  it  is  the  time  average  of  the 
total  mechanical  energy  (kinetic  +  potential)  associated  with  the 
k>h  plant  mode. 

Accordingly  herein,  £*  is  termed  the  plant  mode  energy. 
Similarly,  we  define  the  compensator  mode  energy*: 

A.  =  \  >  (5.4) 

which  may  also  be  interpreted  physically  as  the  electromagnetic 
energy  stored  in  the  inductive/capacitive  elements  of  the  analog 
controller  electronics. 

The  system  dynamics  can  be  understood  in  terms  of  the  energy 
sharing  and  power  flow  between  structure  and  controller.  It  is 
known  quite  generally  that  equations  of  motion  can  be  formed  for 
the  determination  of  the  £*’s  and  £*'s  alone.  To  illustrate  this, 
consider  an  inherently  stable  form  of  the  linear  control  law: 

Kq  —  0,  Kv  = 

(6) 

Fp  —  0,  Fv  ~  k6 

where  k  is  a  real  nonnegative  constant  and  6  €  IR".  This  control 
is  stable  for  all  6  €  IRln,  6  €  IR*  because  iTi+iTi  is  a  Lyapunov 
function  for  the  closed-loop  system.  Applying  the  principles  of  Star 
tistical  Energy  Analysis  |1],  we  obtain  the  following  (approximate) 
equations-of-motion  for  the  plant  and  compensator  modal  energies: 

— Efc  =  -2rtknkEk  +  ]T^(Tfc*(£e  -  Ek)  (7.a) 

e 

k  —  1, .... n 

-  £fc)  (7.6) 

t 

k  =  1, . . . ,  nc 


where  x  =  ^  ^  €  IR2*'  is  the  state  vector  of  a  dynamic  compen¬ 

sator.  Assuming  analog  (continuous  time)  implementation  of  the 
controller,  i  evolves  according  to: 


=  U'klBl  - f,tn' - ; -  (8) 

2  [(r>knk -r  *'{}')*  + (n„  - nt)i) 


o  n  , 

[-/}  -2q/}j 


where  Ko,  Kv,  Fo,  Fv  ve  constant  gain  matrices  and 


f5  =  diag  (A),  A  >0,  »  =  l...,n 
=  diag  (r),},  q,  >  0,  t  =  !...,*» 


«) 


It  is  seen  that  the  generic  linear  compensator  consists  of  a  collection 
of  oscillatory  modes,  as  does  the  plant.  As  will  be  seen,  to  be 
effective,  the  compensator  modal  frequencies,  A.  must  stand  ir  a 
certain  relationship  to  the  plant  frequencies. 

To  probe  some  of  the  limitations  of  linear  compensation  for 
structural  vibration  suppression,  we  first  note  that  the  action  of  the 
compensator  can  be  understood  in  terms  of  its  effect  on  the  energy 
of  vibration.  As  a  measure  of  the  amplitude  of  the  fc'uvibration 
mode,  define: 

Ek  =  ~2  <  >  (5-4 

where  <  >  denotes  a  time  average  over  several  periods  of  vibration. 


(7)  is  a  set  of  power  balance  relations  displaying  the  way  in 
which  the  feedback  gains  meditate  the  exchange  of  energy  among 
the  plant  and  compensator  modes.  (7.a),  for  example,  states  that 
the  rate  of  change  of  the  kth  plant  mode  energy  equals  the  sum  of 
the  power  loss  due  to  dissipation  {-2r}kf7kEk)  and  the  net  power 
flow  from  the  fc* 11  plant  mode  into  all  the  compensator  modes 
(2c  —  £*))•  The  net  power  flows  are  seen  to  be  propor¬ 

tional  to  the  energy  differences  and,  because  of  the  nonnegativity 
of  the  coefficients  <rkt%  power  always  flows  from  the  higher  energy 
mode  to  the  lower  energy  mode.  This  energy  exchange  is  more 
rapid,  the  larger  is  the  power-flow  coefficients  okt.  An  efficient  lin¬ 
ear  controller  design  achieves  its  results  by  making  the  <rkt ’s  as  large 
as  possible,  to  facilitate  energy  transfer  from  plant  to  compensator, 
and  by  choosing  the  17**3  (  the  compensator  modal  damping  ratios) 
somewhat  larger  than  the  rfk's,  thereby  speeding  up  the  dissipation 
of  the  energy  transferred  to  the  compensator. 

Equation  (8)  shows  that  the  power  flow  coefficients  are  in¬ 
herently  nonnegative  and  are  sharply  peaked  functions  of  the  fre¬ 
quency  separation,  1  /?*  -  Ft,  ,  between  plant  and  compensator 
modes.  Thus,  efficient  linear  control  design  (via  Linear-Quadratic- 
Gaussian  design,  for  example)  maximizes  the  <7*/$  by  choosing 
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the  ftk  to  nearly  match  th«  plant  mod*  frequencies.  Thu  feature 
of  quadraiically  optimal  design,  while  it  confers  great  efficiency, 
i*  also  the  source  of  major  limitations.  First,  designed-for  per¬ 
formance  can  be  achieved  only  if  the  plant  modal  frequencies  are 
accurately  estimated  in  advance.  In  any  case,  a  particular  plant 
mod*  exchanges  energy  efficiently  only  with  the  compensator  mode 
that  matches  it*  frequency. 

The  primary  question  addrsssed  her*  is:  Is  it  possible,  by  re¬ 
placing  the  constant  gains  in  (8)  by  functions  of  y  and/or  1,  to 
create  a  nonlinear  compensator  that  achieves  mors  efficient  power 
low  between  plant  modes  and  compensator  modes  -  i.e.,  energy 
exchange  that  is  nearly  independent  of  the  modal  frequency  differ¬ 
ences  and  that  permits  large  power  flow  from  any  one  plant  mode 
simultaneously  to  all  compensator  modes?  In  the  following  Sec¬ 
tions,  w*  propose  a  nonlinear  compensator  design  and  investigate 
the  design  via  both  numerical  simulations  and  analysis.  Although 
the  results  an  by  no  mean*  complete,  these  exploratory  investiga¬ 
tions  indicate  an  affirmative  answer  to  the  above  question. 

W*  find  that,  independently  of  plant  modelling  errors,  the  non¬ 
linear  compensator  provides  very  effective  vibration  suppression. 
Moreover,  the  compensator  can  be  viewed  as  the  interconnection 
of  very  simple  modular  unite  and  its  effectiveness  increases  in  pro¬ 
portion  to  the  number  of  modules.  This  raises  the  question:  Can 
the  proposed  nonlinear  compensator  be  realised  in  a  neural  net? 
Accordingly,  w*  demonstrate,  in  Section  5,  that  the  nonlinear  com¬ 
pensator  can  be  implemented  as  a  neural  net  with  analog  neurons. 

3.  A  Nonlinear  Compensator  for  Structural  Vibration 

Suppression 

With  plant  model  (1)  and  (2),  let  us  consider,  in  place  of  the 
linear  controller  (3),  (4),  (6),  the  nonlinear  controller 


u  -  -<ty*Tij  (9) 


where: 

*T  -  (1, 1 . 1) 

In  effect,  we  have  replaced  *  in  (6)  by  icy  where  y  is  the  sensor 
measurement  signal  (2),  and  *  is  again  a  nonnegative  constant 
whose  magnitude  indicates  the  controller  ‘gain.*  We  now  study 
the  dynamics  of  the  closed-loop  system  defined  by  (1),  (2),  (9), 
and  (10). 

The  intuitive  reasoning  behind  the  choices  (9),  (10)  is  as  fob 
modal  character  of  the  linear 

d 

-2  i\d 

measurement  signal  y  in  order 
to  obtain  a  quadratic  nonlinearity  for  the  compensator  as  a  whole. 
Motivated  by  analogies  with  fluid  dynamic  turbulence,  a  quadratic 
nonlinearity  was  desired  in  order  to  promote  chaotic*  motion  in 
the  closed-loop  system.  This  chaotic  dynamics  endows  the  sig¬ 
nal  y  with  a  smooth,  broad  band  power  spectrum  with  no  spilae* 
or  dominant  harmonics.  The  resulting  broad-band  character  of  the 
‘feedback  gain*,  *y,  is  expected  to  give  rise  to  very  efficient  power- 
flow  from  each  plant  mode  to  all  compensator  modes  in  a  manner 
that  is  largely  insensitive  to  the  precis*  values  of  the  structural 
modal  frequencies- 

To  see  if  the  above  intuitive  motions  were  correct,  we  first  ob¬ 
served  closed-loop  performance  via  ‘brute  force*  numerical  simula¬ 
tion*  for  a  particular  model  of  the  structural  plant.  Specific  results 
and  general  observations  arc  given  in  the  next  Section.  Then  using 
some  of  the  general  empirical  observations  from  the  simulations, 
a  semi-empirical  theory  was  developed  in  the  form  of  a  system  of 

*  In  using  the  term  ‘chaotic  dynamics*  herein,  we  refer  to  the 
operational  definition  given  in  (2|  -  namely,  a  system  exhibits  chaotic 
dynamics  when,  despite  purely  deterministic  initial  conditions  and 
periodic  inputs,  its  measured  response  exhibits  smooth,  continuous 
power  spectra. 


energy  flow  equations  analogous  to  (7).  These  results  are  given  in 

Section  (4). 

3.  Numerical  Simulations  for  an  Example 

For  preliminary  inveatigatioa  of  the  compensator  (9),  (10),  we 
performed  numerical  simulations  for  a  particular  example  of  the 
structural  plant  (1),  (2).  The  example  chosen  is  a  string  extended 
along  x  €  [0,  L\,  held  fixed  at  both  ends  with  uniform  tension  T. 
The  potential  differential  equation  for  the  lateral  deflection,  u>(x,  t) 

U»(0)  w  t »{£)  m  0 

where  p  is  the  constant  lineal  mass  density  and  /(x)  is  the  force 
distribution  due  to  a  single  control  actuator.  The  modal  decompo¬ 
sition  of  this  system  has  the  form: 

u»(*.  t)  =  £  x)wfc(t)  (12) 

* 

D 

*M*)  =  \[jsinkr2 

when,  Assuming  uniform  proportional  damping,  the  modal  coordi¬ 
nates  tv*,  satisfy: 

i 

to*  +  2qf?*th*  +  nlu>k  =  -  dx**(x)/(x )  (13) 

/>  Jo 

Now,  we  nondimensionaliae  variables  so  that  =  1  and  */  = 

and  suppose  that  /(x)  arises  from  a  point  force  actuator 
located  at  x  w  Then: 

t ii*  +  2  lj/7*ti*  +  f7*ui*  =  b*u  (14) 

nk  =  k 

b*  =  sin  kr(  a 

Finally,  assuming  a  collocated  sensor  and  defining  the  plant  state 
as  xT  =  (f?iui|, ...  til, ....  tin),  the  equations  of  motion 

are  found  to  be  identical  in  forms  to  (1)  and  (2)  with: 

O *  =  fc;  k=  l,...,n 

rjk  =  q  =  constant  damping  ratio  (15) 

bT  =  {sin  sin  . . . ,  sin  nxf»] 

With  the  above  expressions  and  for  a  variety  of  choices  of 
d,  H  and  1,  w*  conducted  numerical  simulations  of  the  closed-loop 
systems  consisting  of  (1),  (2),  (9)  and  (10)  with  various  initial 
conditions.  The  qualitative  results  are  not  very  sensitive  to  the 
choice  of  tl  or  the  initial  conditions.  Some  of  these  results  are 
illustrated  in  Figs.  1  and  2,  which  pertain  to  the  case  rjk  =  17  = 
0.002,  da  •  /?*,  *  0.3  and  n  =  20  (so  that  there  are  80  states 

in  the  closed-loop  simulation)  and  with  initial  conditions  such  that 
the  first  mode  has  unit  displacement  and  velocity  and  all  other 
states  are  sero  -i.e.: 

xT(0)  =  (1, 1,0, .  .  .  ,0) 

The  simulations  were  obtained  using  a  fourth-order  Runge-kucta 
integration  routine.  Special  care  has  to  be  exercised  in  selecting  a 
sufficiently  small  integration  time-step,  since,  as  will  be  seen,  y(t) 
exhibit*  very  high  frequency  content  for  sufficiently  large  values  of 

K. 

Fig.  1  shows  time  histones  of  the  displacement  response  of 
the  initially  excited  first  mode  and  the  corresponding  time  histo¬ 
ries  of  the  sensor  measurement  for  three  typical  values  of  x.  For 


lows.  First,  although  we  retain  the 

compensator  -i.e.,  the  term 

are  now  chosen  proportional  to  the 


2,  the  feedback  gains 
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very  small  x,  Che  first  mode  response  shows  a  lightly  damped  peri¬ 
odic  motion  dominated  by  the  first  mode  frequency  and  y(t)  shows 
similar  characteristics.  Slightly  larger  values  of  x  result  in  weakly 
damped  periodic  motions  with  higher  harmonics  of  the  first  mode 
frequency  coming  more  into  play. 

On  the  other  hand  for  x  above  some  critical  threshold  (which 
is  roughly  unity  in  this  example),  the  response  time  histories  ex¬ 
hibit  a  qualitative  change.  As  illustrated  by  the  middle  plots  in 
Fig.  1,  the  initially  excited  mode  drops  dramatically  in  amplitude 
after  a  relatively  brief  period  and  then  is  damped  slowly  thereafter. 
Neither  x,(t)  nor  y(t)  exhibits  any  apparent  periodicities  and  y(t), 
in  particular,  shows  evidence  of  higher  frequency  content.  All  of 
these  tendencies  are  amplified  for  still  larger  values  of  x  (see  the 
bottom  of  Fig.  1). 

Further  insight  into  the  system  dynamics  is  afforded  by  Fig. 
2  which  shows  (under  the  same  conditions  as  in  Fig.  1)  the  time 
histories  of  the  instantaneous  modal  energies  (defined  by  equations 
(5)  but  without  the  time  averaging)  and  the  corresponding  power 
spectra  of  y(t).  For  small  x  (top  part  of  Fig.  2)  energy  sloshes  back 
and  forth  among  the  first  several  plant  modes  and  the  power  spec¬ 
trum  of  y  exhibits  sharp  isolated  spikes.  For  x  above  the  threshold 
(middle  part  of  Fig.  2),  it  is  seen  that  the  rapid  initial  drop  off  of 
the  first  mode  energy  is  accompanied  by  a  redistribution  of  energy 
into  all  the  other  modes,  so  that  after  a  brief  period,  all  the  plant 
and  compensator  energies  are  roughly  equal.  During  the  period 
wherein  modal  energies  are  equalised,  the  total  energy.  By: 

£r  *  (i«i 

*wl  kml 

does  not  appreciably  decline.  E r  is  dissipated  at  a  rather  small 
rate  consistent  with  the  assumed  plant  and  compensator  damping 
ratios  (ij  w  d  »  0.002).  Thus,  the  rapid  decline  in  the  initially  ex¬ 
cited  mode  observed  in  Fig.  1  is  due  not  to  direct  energy  dissipation 
but  to  the  flow  of  the  first  mode  energy  into  all  other  modes.  Ac¬ 
companying  the  modal  energy  equalisation  phenomenon,  the  peaks 
in  the  power  spectrum  of  y  have  broadened  and  coalesced  to  form  a 
continuous  spectrum.  Since  the  spectral  peak  broadening  and  co¬ 
alescence  is  much  larger  than  what  can  be  attributed  to  damping 
and  to  the  finite  time  period  of  the  time  sequence  used  to  calcu¬ 
late  the  spectrum,  it  is  apparent  that  the  system  undergoes  chaotic 
motion  for  x  above  the  critical  value. 

The  above  tendencies  are  strengthened  for  still  larger  values  of 
x  (bottom  part  of  Fig.  2).  Modal  energy  equalisation  occurs  even 
sooner,  the  energy  flow  to  all  modes  occurring  at  nearly  the  same 
rate;  regardless  of  the  relative  values  of  the  modal  frequencies.  The 
power  spectrum  of  y  is  further  smoothed  and  broadened.  Indeed, 
the  spectrum  of  y  is  nearly  constant  over  the  whole  frequency  band 
occupied  by  the  plant  and  compensator  modal  frequencies. 

The  above  findings  tend  to  confirm  the  heuristic  insights  used 
in  constructing  the  design  (9),  (10).  In  particular  the  quadratic 
nonlinearities  in  the  compensator  do  trip  the  system  into  chaotic 
motion,  resulting  in  a  broad-band  spectrum  for  y(t)  and  very  ef¬ 
ficient  energy  flow  among  all  the  plant  and  compensator  modes, 
even  for  modes  having  widely  separated  frequencies. 

Surveying  all  the  simulation  results,  we  obtain  additional  ob- 
servatione  regarding  time-averaged  correlation,  autocorrelation 
functions  and  power  spectra  that  prove  useful  in  constructing  a 
semi-empirical  theory  of  the  dynamics  of  the  nonlinear  compen¬ 
sator.  With  regard  to  correlations  and  autocorrelations,  we  have: 

0.1  For  t  larger  that  ~10  lowest  mode  periods,  the  separate 
modal  coordinates  zia.xjs;  k  —  1 . n  are  approxi¬ 

mately  uncorrelated  (in  the  sense  of  time  averaging) 


In  addition,  we  have  observations  concerning  the  changing 
character  of  the  power  spectra  of  the  plant  modal  velocities,  ij* 
*  =  1, . . . ,  n,  as  x  increases.  We  denote  the  power  spectrum  of 
X2S  by  s»„(w).  Noting  that  the  correlation  coefficient,  p„,.(r),  is 
the  inverse  Fourier  transform  of  5„,(w)/ /“  dwS,„,  the  following 
observations  also  have  direct  import  for  correlation  functions: 


0.2  For  very  small  x,  5,u  exhibits  isolated  spikes  at  the  fre¬ 
quencies  of  excited  modes  having  half-power  widths  equal 
to  2/}kn 

0.4  For  larger  x,  near  the  critical  threshold  value  x., 

shows  spikes  at  many  additional  modal  frequencies.  The 
width  of  the  spikes  grows  to  —  x£^  [Ey  given  by  (16)J. 
The  value  of  xc  appears  to  be  roughly  -Ij-iw,  where  6u 

bt 

is  the  minimum  separation  between  plant  modal  frequen¬ 
cies. 


0.5  For  x  >>  x„,  the  spikes  in  S**k  coalesce  into  a  smooth, 
broad-band  spectrum,  which  is  approximately  constant 
up  to  some  frequency  /?*  and  then  drops  off  rapidly  at 
higher  frequencies.  The  value  of  17*  is  roughly  *Ey  . 

4-  A  Semi-Empirical  Theory  for  the  Energy  Dynamics  of 
the  Nonlinear  Compensator 


Here  we  use  the  simulation  results  and  corresponding  observa¬ 
tions  discussed  in  the  last  Section  to  construct  a  semi-empirical  the¬ 
ory  for  the  nonlinear  compensator  (9),  (10).  The  theory  takes  the 
forma  of  approximate  equations-of-motion  for  the  time-averaged 
modal  energies  analogous  to  (7).  These  equations  are  then  used 
to  deduce  various  qualitative  phenomena  and  provide  a  few  useful 
design  guidelines. 

Space  limitations  preclude  the  full  derivations,  which  will  be 
given  elsewhere.  Here  we  attempt  merely  to  sketch  the  develop¬ 
ment. 


The  first  step  is  to  form  equations  of  motion  for  the  ‘second 
moment  matrix*,  Q  =  ^ j  (xT,  iT),  of  the  full  closed-loop  sys¬ 
tem  (1),  (2),  (9)  and  (10).  We  then  manipulate  the  equations  to 
eliminate  the  cross-correlation  terms  in  favor  of  the  modal  mean- 
squares  xj,  xj,;  k  =  i...n,  m  =  l...nc,  apply  the  time  av¬ 
eraging  operator  and  employ  a  perturbation  expansion  approach 
to  obtain  equations  approximately  valid  for  small  x.  Neglecting 
terms  of  order  xJ  or  smaller,  we  obtain  the  following  equations 
for  the  time-averaged  plant-mode  energies,  £*,  and  compensator 
mode  energies  £*: 


Ek  *  -2q*/7*£*  +  ^d*r(£f  -  £*),  k=l...n 

t 

&k  —  —  2^*/5»£*  +  ^  &ik(Et  —  £*),  k  =  l...ne 

e 

where: 

dkt  =  x3»2  f  dr  <  y(t)y(t  -  r)  >  7  «(r) 

Jo 

7»t(0  =  ^*"i’,‘°*'r’’'rt‘|,|cos(tv*  -w()t  +  cos(w*«()t) 
and  where  w*,w*  denote  the  damped  natural  frequencies: 


(18) 


(19) 


“t  =  ftk\J  1  -  nl 

u>*  =  -  7* 


(20) 


0.2  Again  for  t~  to  lowest  mode  periods,  the  autocorrelation 
coefficient  of  y: 

p*(t.  *)  =  <  y(t)y(t  -<■)>/<  v2(<)  >  (i?) 

is  approximately  independent  of  t  (i.e.,  it  is  weak-sense 
stationary) 


Although  developed  for  small  k,  equations  (Id)  appear  to  give 
correct  results  even  for  large  *.  This  may  be  due  to  the  semi* 
er  pineal  manner  in  which  explicit  expressions  for  the  coefficients 
dkt  are  derived,  as  discussed  in  the  following. 

It  remains  to  express  dkt  explicitly  in  terms  of  the  modal  ener¬ 
gies.  To  do  tjiis  we  use  the  empirical  observations  0.1  through  0.5 
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given  in  the  previous  section.  Pint,  tad  immsdUt*  consequence  of 
0.1  tad  0.2  it: 

<  »(0v(‘  -»•)>»  53  k*  <  *l*(0  >  P..(r)  (21) 

k 

when 

*•.(*■)  ■<*»(<)*»(< -*)  >/<  *i*M  >  (22) 

However,  neglecting  higher-order  ttrmt  is  a;  <  *j*  >  at  £»(<)• 
Using  this  approximation  ia  (21)  tad  substituting  tht  malt  iato 


*u  as  *J4*  53  (23.a) 

«* 

f  <tr/>Mm(rhu(r)  (23.4) 

Jo 


Next  we  deduce  the  form  of  by  coaiideriug  two  limitirg 
cteet:  very  email  <c  tad  very  large  *.  Ia  the  cate  of  very  tmall  x, 
0.3  impliet: 

P._(r)  «  •'""“'cot  wmr 

Thit  caa  be  rabatituted  into  (23.b)  aad  /Vet  evaluated  directly. 
For  large  it,  we  ute  0.4  and  0.5  together  with  dimeutional  aaaiyeit 
to  deduce  the  atymptotic  form  of  rs>a(.  At  the  latt  etep,  an  expree- 
tion  for  /*,*«  it  devittd  which  correctly  reduce!  to  the  exprettiont 
derived  for  the  two  limiting  caret.  The  final  retult  it: 


(24.a) 


/m*£  *  j(4m  +  Ik Ok  +  fiA) 

{ - i _ 

(An  +  +  fiA)3  +  («n  +  W*  -  Wj)1 

_  1 _ 

(^m  +  +  fit  fit)3  +  (u,.  -U»  +  WtP 

_ 1 _ 

(An  +  +  fit  f)t)J  +  (ua  +  +  <5t)J 

_ 1 _ 

(An  +  IkOk  +  fit  l)t)a  +  (wm  -  W*  -  Wt)3 
4m  —  *?ffm  +  kE^ 


(24.4) 

(2d-*) 


where  wa,wa  and  £T  are  defined  in  (20)  and  (16),  respectively. 

Equations  (18)  and  (24)  constitute  a  doted  lystem  of  equa¬ 
tions  approximately  describing  the  dynamics  of  the  modal  energies 
of  the  plant  and  compensator  and  may  now  be  used  to  deduce 
various  properties. 


First,  it  should  be  noted  that  (18)  are  of  the  tame  structure 
at  (7).  The  power  Bow  from  any  one  mode  to  all  other  modes 
it  again  proportional  to  the  energy  differences  and  because  the 
coupling  coefficients  dat  are  all  intrinsically  positive,  poster  always 
flows  from  the  more  energetic  to  the  lets  energetic  mode.  These 
features  ensure  that  the  system  will  be  driven  toward  equalisation 
of  modal  energies  with  the  time  scale  for  equalisation  being  dictated 
by  the  magnitude  of  kE£.  This  is  consistent  with  the  qualative 
observations  of  the  latt  section.  Furthermore,  (24)  shows  that  for 


sufficiently  large  *;  rmat  approaches  the  uniform  limit 
that: 


kE. 


c 


<  y3 


> 


Thus,  in  contrast  to  (8),  there  is  strong  power  flow  from  any  one 
plant  mode  to  all  other  compensator  modes,  regardless  of  the  rel¬ 
ative  values  of  the  plant  compensator  frequencies.  This  efficient 
energy  sharing  is  a  consequence  of  the  nonlinearity  introduced  in 
the  compensator  design  (9)  and  (10). 

As  a  last  topic,  we  use  (18)  and  (24)  to  obtain  simple  quantita¬ 
tive  estimates  of  the  speed  with  which  structural  vibration  energy 


can  be  drained  away  to  the  compensator  via  nonlinear  design.  For 
thit  purpose,  consider  the  cate  wherein  a  set  of  n<  structural  modes 
are  directly  excited  and  it  it  desired  to  reduce  the  vibration  energy 
of  these  modes  because  of  their  determinout  impact  on  system  per¬ 
formance.  We  estimate  the  doted- loop  response  by  taking  account 
of  the  interactions  between  the  compensator  modes  and  only  these 
directly  excited  modes.*  For  simplicity,  suppose  that  all  nd  modes 
are  initially  excited  to  the  tame  energy: 

£a(0)  «  Eo  Vfc  =  l...m 

aad  that  the  compensator  is  initially  quiescent;  Le.,  £jt0>  »  0  V*. 
Also  let  us  estimate  the  magnitude  of  all  the  modal  influence  co¬ 
efficients  by  some  average  value  {,  Le.,  4j  at  l3  Vfc.  Finally,  since 
flexible  structure  vibration  control  is  our  motivating  application, 
see  assume  tmall  damping  ratios  for  both  plant  and  compensator; 
i.e.,  to,  «  1.0,  fia  <<  1.0. 

With  the  above  conditions,  supp~>ee  that  [0,412]  is  the  fre¬ 
quency  band  encompassed  by  all  the  initially  excited  modes.  We 
need  only  choose  the  n,  compensator  frequencies,  /)»;  k  =  1 .  . .  n„, 
somewhere  in  this  band,  because,  as  (24. b)  shows,  the  choice  of  x 
such  that 

x£J  >  34/1  (25) 

i mres  that  all  the  /'mac’s  reduce  to  the  same  value,  -4,  inde- 

pendently  of  the  compensator  frequencies.  The  design  choice  (25) 
implies,  by  use  of  (24.a),  that: 


°kt  — 


(26) 


Using  this  result  and  the  fact  that  during  the  initial  period  of  en- 
ergy  equalisation,  Em  can  be  estimated  by  £T,  equations  (18) 
become: 


Eh  =  -2 vkflkEk  +  kE*P  £(£c  -Ek);  k  =  1 . . .  n„ 

e«i 

£  —  —  2>7k/\£fc  +  —  ■£*)»  k  as  1 . . .  ne 


Ek(Q)  =  E0,  £k(0)  =  0  V* 

Note  that  the  above  equations  imply: 


—  =  -2  VkftkEk  + 


(28) 


JEt(O)  =  ndE o 

Evidently,  the  total  energy  is  dissipated  over  a  time  scale  of  or¬ 
der  ^  or  ^  (for  some  k ).  By  virtue  of  our  small  damping  ratio 
assumption  and  design  choice  (25),  this  time  scale  is  much  longer 
than  the  time  scale  over  which  the  initial  energy  redistribution 
takes  place.  For  investigation  of  the  initial  period  of  energy  equal¬ 
isation,  therefore,  we  may  treat  Et  as  a  constant  2:  n.iEo  and 
neglect  the  damping  terms  involving  rjk  and  in  (27).  With  these 
approximations,  one  immediately  obtains  the  following  solutions 
for  the  relatively  brief  initiai  time  period  over  which  energy  redis¬ 
tribution  occurs: 

Ek  =  Eo - 2_(1  -  e',/T*);  k  =  1 ...  n, 1 

nEnd  +  n'  (29) 

£t  =  -2i£_(I_e-/T  ,)i  k=  1  n 

■7’ 


where: 

Tk  =  -  T  - -  (301 

*i'2Efnti{n,j  n.) 

•  This  simplification  actually  results  in  an  overestiniauon  of  the 
energies  resident  in  the  nj  excited  modes. 
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Thus,  afar  a  tim*  period  of  order  T„  the  energy  initially  rc*  id- 
mi  in  the  plant  structural  vibration  ia  drained  away  to  the  compen- 
eator  modes  and  all  modal  energies  are  approximately  equalised 
to: 

f.  _  &T 


t>T„:  Ekn 


+  n. 


,  4s» 


n S  +  ". 


VJfc 


(31) 


Now  for  t  >>  Tm,  we  can  characterise  the  evolution  of  the 
total  energy,  By,  by  using  (31)  in  (28)  to  obtain: 


t  » 


r„:  ^By  *  -2  n»nk  +  £  Er 


(32) 


Thus,  over  large  time  scales,  the  total  energy  is  damped  exponen¬ 
tially  with  an  equivalent  damping  ratio  that  it  a  weighted  average 
of  ieli  the  plant  and  compensator  damping  ratios. 

Another  noteworthy  feature  of  the  above  results  is  that  both 
the  time  period,  Tm  (equation  (30)),  needed  (or  equalisation  of  en¬ 
ergies  and  the  value  to  which  the  modal  energies  are  equalised 
(equations  (31))  are  univeraely  proportional  to  the  number  of  com¬ 
pensator  modes.  Thus,  especially  effective  vibration  suppresaion 
can  be  achieved  by  the  nonlinear  compensator  if  the  number  of  its 
states  can  be  made  very  large.  This,  suggests  the  question:  Does 
the  compensator  (9),  (10)  have  a  simple  repetitive  structure  can 
this  structure  be  implemented  as  a  neural  net  containing  a  large 
number  of  analog  neurons? 

5.  A  Neural  Not  Realisation  of  the  Nonlinear 
Compensator 

Equations  (9)  and  (10)  may  be  rewritten  to  reveal  that  the 
control  signal  u  is  the  sum  of  nc  components: 


u  =  £* 


(33) 


is  a  Lyapunov  function  for  the  systems  (1),  (2),  (33),  (38),  since: 

dJ  *  _  *• 

»  -2  ^  ’Tfartfc*2fa  -  53  faAk(g{tk)Sk  +  ?(&fa)$fa)  (37) 

fa-1  fa-1 

Thus  the  controller  (33),  (35)  it  inherently  stable. 

Summary  end  Conclusion 

in  this  paper,  we  explored  a  novel  type  of  nonlinear  dynamic 
controller  design  that  was  originally  motivated  by  certain  issues 
in  the  problem  of  vibration  suppression  in  flexible  structures,  fa 
cause  of  the  particular  form  of  quadratic  noshaearitiee,  the  con¬ 
troller  provides  extremely  efficient  energy  exchange  mechanisms 
capable  of  rapidly  draining  energy  sway  from  the  structural  plant, 
la  addition  the  controller  was  shown  to  consist  at  the  interconnec¬ 
tion  of  simple  repetitive  modules  and  its  effectively  increases  with 
the  number  of  such  modules  and  its  effectively  increases  with  the 
number  of  such  modules.  These  features  motivate  the  implementa¬ 
tion  of  the  nonlinear  controller  via  the  neural  network  architecture 
explored  in  the  last  Section.  With  the  advent  of  appreciate  ana- 
log  neural  net  hardware,  this  suggested  nonlinear  compenaation 
scheme  could  offer  a  very  effective  means  of  vibration  suppression. 
For  example  suppoee  tome  10  structural  modes  are  significantly 
excited  and  must  be  suppressed  and  we  employ  the  neural  net 
controller  involving  a  modest  number  of  neurons,  say  2000.  Then 
*  10,  e,  »  10s  and  it  is  seen  from  (31)  that  the  vibrational 


of  ite 


where  each  u*  is  the  output  of  a  simple  nonlinear  oscillator: 


energy  of  each  excited  mode  is  quickly  reduced  to 
.  .  .  si  +  n, 

initial  value  -  a  reduction  of  more  than  a  hundredfold. 
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(34) 


Thu*,  the  nonlinear  compensator  is  a  combination  of  simple  repet¬ 
itive  modules  and  this  suggests  that  it  can  be  efficiently  imple¬ 
mented  as  a  neural  network. 

To  see  that  this  is  the  case,  define  y(  )  to  be  some  antisym¬ 
metric  sigmoidal  function  such  that  g' (x)  is  maximum  at  *  *  0, 
where  $'(0)  =  1.  Then  consider  the  replacement  of  (34)  by: 

Ufa  =  -Kyg(vh) 

tug) -«•(:) 

(35) 

The  above  equations  essentially  reduce  to  (34)  for  small  signal 
amplitudes  such  that  j(x)  as  x. 

It  is  now  easy  to  see  that  (33),  (35)  are  equivalent  to  a  system 
formed  by  interconnection  of  analog  neurons,  of  the  form  given 
by  Hopfield,  illustrated  here  in  Fig.  3.  Using  such  neurons,  we 
form  neuron  pairs  in  the  manner  shown  in  Fig.  4. a,  such  that 
each  neuron  pair  implements  one  compensator  mode.  We  finally 
interconnect  n„  neuron  pairs  as  shown  in  Fig.  4.b  to  obtain  a 
system  completely  equivalent  to  (33),  (35). 

It  is  seen  that  this  as  an  implementation  involving  very  sparse 
neuronal  interconnections.  Note  also,  that  just  as  ^  £*(*?*  + 
*»)  +  j  +  *»)  *  a  Lyapunov  function  for  the  closed  loop 

system  (1),  (2),  (9),  (10);  so  too,  the  quantity 


■  "  "■  /  /-<.(« I  /■*«!*  1  \ 

J=2  + £■'(/>  gMdx  +  J0 


Fig.  1:  Time  histories  of  the  initially  excited  modal  displacement 
and  the  sensor  measurement  for  fl  —  fit,,  9s  =  0.002  and 
various  values  of  k. 


(36) 
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Fig.  2:  Tim*  histories  (for  th*  urn*  condition!  u  Fig.  1)  of  tht 
instantaneous  model  energies  (£„  =  ^(ij„  +  i^),  £«  • 
3(if„  +  i|.))  and  the  corresponding  power  spectra  of  the 
tensor  measurement. 


f,  .  bund  i«f«i  im*  w  a»e«sa  * 

Fig.  3:  Basic  structure  of  an  analog  neuron  following  Hopfield 

(3). 


<b) 


Fig.  4  Implementation  of  controller  (33),  (3$)  via  a  neural  net¬ 
work: 

(a)  fundamental  neuron  pair,  (b)  overall  architecture. 
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Abstract 

Real  parameter  uncertainly  and  phase  information  play  a  hey  role  in 
the  analysis  and  synthesis  of  robust  controllers  for  lightly  damped  flexible 
structures.  In  this  paper  we  discuss  the  ramifications  of  this  issue  as  it  af¬ 
fects  achievable  performance  in  structural  control.  In  this  regard  we  review 
the  state  of  knowledge  in  addressing  real  parameter  and  phase  wanes.  The 
discussion  is  illustrated  by  examining  robust  controllers  designed  for  the 
ACES  structure  at  Marshall  Space  Flight  Center.  These  controllers  were 
designed  by  means  of  the  Maximum  Entropy  generalised  LQG  methodol¬ 
ogy 

1.  Introduction 

Traditionally,  spacecraft  control-system  designers  have  been  primar¬ 
ily  concerned  with  controlling  rigid  body  attitude  modes,  while  avoiding 
the  excitation  of  flexible  body  dynamics.  As  performance  requirements 
become  more  stringent  and  spacecraft  become  larger,  control- system  de¬ 
sign  must  explicitly  encompass  flexible  dynamic  modes  so  as  to  actively 
suppress  undesired  structural  vibration.  Furthermore,  for  complex  space¬ 
craft,  multi-input  multi-output  controllers  with  significant  bandwidth  may 
be  required. 

Since  structural  modeling  and  identification  of  large  flexible  structures 
in  a  1-g  environment  possess  inherent  limitations,  one  of  the  key  issues  in 
structural  control  is  robustness.  Although  robust  control  has  undergone 
intensive  development  in  the  past  two  decades,  there  remain  aspects  of  ro¬ 
bust  control  that  are  relevant  to  structural  control  and  that  are  largely 
unresolved.  These  aspects  are  the  role  of  real  parameter  uncertainty  and 
phase  information.  The  purpose  of  this  paper  is  to  examine  the  impact  of 
these  issues  on  structural  control,  their  interrelationship,  and  their  mani¬ 
festation  within  the  analysis  and  synthesis  of  feedback  systems. 

2.  Phase  Stabilisation  Versus  Gain  Stabilisation 

From  a  classical  control-design  point  of  view,  the  issues  of  real  param¬ 
eter  uncertainty  and  phase  information  are  manifested  in  the  fundamental 
concepts  of  gain  and  phase  stabilisation.  In  terms  of  gain  stabilisation, 
stability  of  a  single-input  single-output  closed-bop  system  is  insured  by 
designing  the  controller  so  that  the  magnitude  of  the  loop  transfer  func¬ 
tion  is  less  than  unity  in  frequency  regimes  in  which  the  phase  is  either 
known  to  be  near  180°  or  is  highly  uncertain.  In  terms  of  phase  stabilisa¬ 
tion,  stability  is  achieved  by  insuring  that  the  phase  of  the  loop  transfer 
function  is  well  behaved  where  the  loop  transfer  function  has  gain  greater 
than  unity.  Roughly  speaking,  phase  stabilization  can  be  used  to  allow 
high  loop  gain  and  thus  achieve  high  performance  in  frequency  regimes  in 
which  sufficient  phase  information  is  available,  whereas  gain  stabilisation 
(e.g.,  rolloff)  is  needed  to  insure  stability  where  the  phase  of  a  system  is 
very  poorly  known.  For  further  discussion  of  the  distinction  between  phase 
and  gain  stabilization,  see  (1). 

3.  Structured  Real  Parameter  Uncertainty  Versus  Unstructured 
Complex  Parameter  Uncertainty 

A  variety  of  approaches  have  been  proposed  for  addressing  uncertainty 
in  the  synthesis  of  robust  controllers.  These  include  Hw  synthesis  [2- 
$],  quadratic  Lyapunov  functions  (6-9),  and  the  structured  singular  value 
|l0,il|.  All  of  these  methods  effectively  treat  the  uncertain  parameters  as 
complex  quantities,  and  are  thus  conservative  with  respect  to  real  param¬ 
eter  uncertainty.  If  the  uncertain  parameters  are  known  to  be  real,  then 
special  techniques  are  required  to  avoid  conservatism  |12-19j. 

To  illustrate  the  conservatism  of  H„  theory  in  the  presence  of  phase 
information,  it  need  only  be  noted  that  |e>*;  •>  1  regardless  of  the  phase 
angle  4-  Indeed,  any  robustness  theory  based  upon  norm  bounds  will  suf¬ 
fer  from  the  same  shortcoming.  Of  course,  svery  real  parameter  can  be 
viewed  as  a  complex  parameter  with  phase  4  «  CP  or  4  ■  1*0°.  Since  the 
existence  of  a  single  Lyapunov  function  for  a  norm-bounded  uncertainty 
class  is  equivalent  to  a  small-gain  condition  !9j,  much  of  Lyapunov  theory 
exhibits  a  similar  conservatism. 

In  structural  modeling  via  finite  element  models,  uncertainty  in  the 
mass,  damping,  and  stiffness  matrices  is  unavoidable.  If  the  man  and  stiff- 
ness  matrix  uncertainty  is  modeled  as  complex,  unstructured  perturbations, 
then  the  damping  matrix  is  effectively  perturbed  as  well.  Indeed,  damping 
is  sometimes  modeled  as  a  complex  stiffness  20,  p.  194).  Difficulty  vises 
when  stiffness  uncertainty  is  large  relative  to  damping  uncertainty,  in  which 
case  complex  stiffness  uncertainty  corresponds  to  a  physically  unrealizable 
unstable  plant  model. 

4.  Phase  Information  and  Positive  Ileal  Transfer  Functions 

Phase  information  plays  a  fundamental  role  in  structural  control.  For 

illustration,  consider  a  flexible  structure  with  a  colocated  rale  sensor/force 
actuator  pair  and  assume  these  devices  arc  ideal.  For  such  a  system  the 
transfer  function  from  the  actuator  to  sensor  is  known  to  be  positive  real, 
that  is,  to  have  phase  lying  between  90"  and  -90"  '21,22!.  In  a  negative 
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feedback  configuration,  a  controller  for  this  plant  that  is  strictly  positive 
real  cannot  destabilise  the  system  since  the  loop  transfer  function  has  phase 
leas  than  - 180°  over  all  frequencies.  Hence  such  a  control  system  will  be  un¬ 
conditionally  robust  to  uncertainties  in  both  natural  frequencies  and  damp¬ 
ing.  Of  course,  these  observations  assume  perfect  sensors  and  actuators  so 
as  not  to  introduce  additional  phase  lag.  If  the  sensors  and  actuators  do 
have  significant  dynamics,  then  the  feedback  law  must  be  chosen  so  that 
the  transfer  function  consisting  of  the  cascaded  sensor,  compensator,  and 
actuator  dynamics  is  strictly  positive  real.  If,  in  practice,  positive  madness 
can  only  be  enforced  over  a  limited  frequency  band,  then  loop  gain  rolloff  is 
required  when  phase  lags  or  phase  uncertainties  reach  unacceptable  levels. 

By  exploiting  the  stability  guarantee  due  to  the  interconnection  of  pos¬ 
itive  real  MIMO  systems,  robust  positive  real  controllers  have  been  studied 
for  structural  control  [23-31).  A  related  approach  involves  using  H«,  de¬ 
sign  in  conjunction  with  the  bilinear  transformation  (32).  By  using  a  Riccati 
equation  to  enforce  a  positive  real  constraint,  robust  controllers  for  positive 
real  uncertainty  were  obtained  in  (33).  Related  results  appear  in  (34|. 

Alternative  approaches  to  including  phase  information  in  analysis  and 
synthesis  include  [35-38).  Ref.  [39)  extends  the  gain  envelope  approach  of 
(40,4l[  to  include  a  phase  envelope  as  welL  These  envelopes  are  charac¬ 
terised  by  real  parameters  whose  effect  can  then  be  addressed  using  real 
parameter  robustness  techniques. 

An  alternative  approach  to  exploiting  phase  information  is  based  on 
the  concept  of  structured  covariance  matrix.  Roughly  speaking,  robustness 
is  not  guaranteed  by  means  of  a  Lyapunov  function  or  covariance  bound  [8|, 
but  rather  by  means  of  a  covariance  matrix  whose  structure  is  insensitive 
to  a  given  class  of  plant  perturbations.  This  concept  provides  the  basis 
for  the  generalised  LQG  synthesis  technique  known  as  Maximum  Entropy 
design  (42-48). 

5.  An  Illustrative  Example  Using  Maximum  Entropy  Synthesis 

The  ACES  experimental  testbed  is  located  at  NASA  Marshall  Space 
Flight  Center.  The  basic  teat  article,  a  spare  Voyager  Astromast,  is  a 
deployable,  lightweight  (about  5  pounds),  lightly  damped  beam,  approxi¬ 
mately  45  feet  in  length.  The  Astromast  is  symmetric  with  a  triangular 
cross  section.  Three  longerons  form  the  converse  of  the  beam  and  extend 
continuously  along  its  full  length.  The  cross  members,  which  give  the  beam 
its  shape,  divide  the  beam  into  91  sections  each  having  equal  length  and 
mass  and  similar  elastic  properties.  When  fully  deployed,  the  Astromast 
exhibits  a  longitudinal  twist  of  approximately  260  degrees. 

The  ACES  configuration  consists  of  an  antenna  and  counterweight  legs 
appended  to  the  Astromast  tip  and  the  pointing  gimbal  arms  at  the  As¬ 
tromast  base.  The  addition  of  structural  appendages  creates  the  ‘nested* 
modal  frequencies  characteristic  of  large  space  structures.  Overall,  the 
structure  is  very  flexible  and  lightly  damped.  It  contains  many  closely 
spaced,  low  frequency  modes  (more  than  40  modes  under  10  Hz).  The 
ACES  configuration  is  dynamically  traceable  to  future  space  systems  and 
is  particularly  responsive  to  the  study  of  line-of-sight  (LOS)  issues. 

The  goal  of  the  control  design  is  to  position  the  laser  beam  in  the  cen¬ 
ter  of  the  detector.  The  disturbances  were  chosen  to  be  position  commands 
to  the  Base  Excitation  Table  (BET).  The  BET  motion  is  regulated  by  an 
analog  controller  which  allows  any  type  of  BET  movement  within  the  fre¬ 
quency  limitations  of  the  hydraulic  actuation  system.  In  the  discussion  that 
follows  we  will  consider  only  one  single-input,  single-output  loop  involving 
AGS-X,  the  x-torque  of  the  Advanced  Gimbal  System,  and  BGYRO-X,  the 
rotational  rate  of  the  base  gyro. 

For  the  AGS-X  to  BGYRO-X  loop  a  model  was  developed  by  using  the 
Eigensystem  Realisation  Algorithm.  The  ERA  model  was  compared  with 
the  frequency  response  functions  (FRF's)  derived  from  the  test  data.  The 
ERA  model  matched  the  FRF  data  fairly  closely  in  msgnitnde  although 
the  modal  frequencies  do  not  exactly  coincide.  The  ERA  model  differed 
even  more  from  the  FRF’s  in  phase. 

Control  design  for  LOS  performance,  initially  performed  by  using  stan¬ 
dard  LQG  techniques,  required  penalising  only  the  modes  less  than  3  Hz. 
Thus,  high  performance  controllers  were  limited  to  having  gain  only  at  the 
inodes  less  than  3  Hz.  To  avoid  destabilizing  the  two  higher  frequency 
modes  of  the  ERA  model,  the  LQG  controllers  contained  notches  at  the 
two  corresponding  frequencies. 

The  LQG  controllers  tended  to  be  very  sensitive  to  the  phase  uncer¬ 
tainty  in  the  performance  region,  the  frequency  interval  from  DC  to  3  Hz. 
They  also  were  very  sensitive  to  the  frequency  uncertainty  in  the  two  higher 
frequency  modes.  This  control  problem  thus  provides  an  excellent  real-life 
example  of  phase  uncertainly  and  real  parameter  (in  this  case  frequency) 
uncertainty. 

Robust  control  design  was  performed  using  the  Maximum  Entropy 
(ME)  approach  '49]  This  approach  allows  the  designer  to  directly  account 
for  real  parameter  uncertainty  142-481.  Figure  1  describes  the  influence  of 
ME  uncertainty  design  on  the  phase  of  a  full-order  compensator  in  the  per¬ 
formance  region.  The  phase  of  the  LQG  compensator  varies  widely  (and 
wildly)  over  this  frequency  interval,  implying  that  the  Nyquist  plot  of  the 
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coma  ponding  bop  transfer  function  encircles  the  origin  several  timaa.  Aa 
oat  «nU  expect,  these  daaigu  war*  nonrobust  and,  in  fact,  wan  nnatabla 
whan  implemented.  Howavar,  tba  ME  daaigna  bacama  poaitiva  real  in  tba 
parformanca  ngion  landing  toward  rata  feedback.  Tbua,  tha  ME  daaigna 
provided  tba  seeded  stability  robaatnaaa  in  tba  parformanca  region.  In  ad¬ 
dition,  tba  ME  daaigna  roboatifiad  tha  LQG  controllar  notchaa  by  inenaaing 
both  tba  width  and  depth  of  tha  notchaa  (Figun  2). 
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Abstract 

In  many  applications  of  feedback  control,  phase  information  is  available  concerning  the  plant 
uncertainty.  For  example,  lightly  damped  flexible  structures  with  colocated  rate  sensors  and  force 
actuators  give  rise  to  positive  real  transfer  functions.  Closed-loop  stability  is  thus  guaranteed  by 
means  of  negative  feedback  with  strictly  positive  real  compensators.  In  this  paper,  the  properties 
of  positive  real  transfer  functions  are  used  to  guarantee  robust  stability  in  the  presence  of  positive 
real  (but  otherwise  unknown)  plant  uncertainty.  These  results  are  then  used  for  controller  synthesis 
to  address  the  problem  of  robust  stabilization  in  the  presence  of  positive  real  uncertainty.  One  of 
the  principal  motivations  for  these  results  is  to  utilize  phase  information  in  guaranteeing  robust 
stability.  In  this  sense  these  results  go  beyond  the  usual  limitations  of  the  small  gain  theorem  and 
quadratic  Lyapunov  functions  which  may  be  conservative  when  phase  information  is  available.  The 
results  of  the  paper  are  based  upon  a  Riccati  equation  formulation  of  the  positive  real  lemma  and 
thus  are  in  the  spirit  of  recent  Riccati-based  approaches  to  bounded  real  (Hoc)  control. 
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1.  Introduction 


In  many  applications  of  feedback  control,  phase  information  is  available  concerning  the  plant 
uncertainty.  For  example,  lightly  damped  flexible  structures  with  colocated  rate  sensors  and  force 
actuators  give  rise  to  positive  real  transfer  functions.  Closed-loop  stability  is  thus  guaranteed  by 
means  of  negative  feedback  with  strictly  positive  real  compensators.  This  principle  has  been  widely 
used  to  design  robust  controllers  for  flexible  structures  [1-10]. 

The  salient  features  of  positive  real  transfer  functions  is  that  they  are  dissipative  and  phase 
bounded  [11-25].  Hence  the  feedback  interconnection  of  positive  real  transfer  functions  is  guaran¬ 
teed  to  be  stable  without  requiring  that  a  small  gain  condition  be  satisfied.  Positive  real  design  is 
thus  potentially  less  conservative  than  bounded  real  (Hoo)  design  in  the  presence  of  phase  informa¬ 
tion. 

In  this  paper  we  utilize  properties  of  positive  real  transfer  functions  to  develop  new  conditions 
for  robust  stability  and  robust  stabilizability.  Although  related  results  have  been  developed  previ¬ 
ously  [26-30],  this  paper  goes  beyond  earlier  work  by  exploiting  a  Riccati  equation  formulation  in 
the  spirit  of  recent  advances  in  Hoo  synthesis  [31-37],  This  is  done  in  two  different,  but  equivalent, 
ways.  First  we  show  that  the  Riccati  equation  used  to  enforce  an  Hoo  constraint  can  be  transformed 
to  yield  a  different  Riccati  equation  that  enforces  a  positive  real  constraint  (Theorem  3.2).  Alter¬ 
natively,  we  show  that  the  same  Riccati  equation  can  be  obtained  by  manipulating  the  conditions 
of  the  positive  real  lemma  (Proposition  3.3).  Many  of  the  techniques  and  transformations  used  in 
these  steps  are  due  to  [17],  which  contains  an  extensive  treatment  of  positive  real  and  bounded  real 
transfer  functions. 

Once  the  Riccati  equation  that  enforces  positive  realness  has  been  derived,  robust  stability  can 
be  guaranteed  for  a  class  of  perturbations  involving  an  arbitrary  constant  positive  real  matrix  (see 
the  set  U  defined  by  (4.6)  and  Theorem  4.1).  The  modeling  of  matrix  uncertainty  by  means  of  a 
"fictitious”  feedback  loop  (linear  fractional  transformation)  is  directly  analogous  to  the  small  gain 
(Hoo)  parameter  uncertainty  model  of  [37].  In  our  case,  however,  the  class  of  uncertainties  includes 
a  phase  constraint  rather  than  a  small  gain  condition  (see  Remark  4.1). 

Having  enforced  robust  stability  for  positive  real  uncertainty,  we  then  proceed  in  Section  5 
to  give  necessary  and  sufficient  conditions  for  robust  stabilizability  in  terms  of  a  pair  of  coupled 
algebraic  Riccati  equations  (Theorem  5.1).  A  robustly  stabilizing  feedback  gain  is  then  given  in 
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terms  of  the  solutions  to  the  Riccati  equations.  The  stabilizability  result  is  first  stated  for  static 
output  feedback  and  then  specialized  to  the  case  of  full-state  feedback. 

Finally,  we  close  the  paper  by  discussing  connections  between  the  positive  real  uncertainty 
modeling  approach  of  this  paper  and  the  Maximum  Entropy  approach  to  robust  control  design  of 
[10,38-43]. 

Notation: 

IR,IRrx*  real  numbers,  r  x  s  real  matrices 

Ir ,  /;  (  )T,  (  )*  r  x  r  identity  matrix;  transpose,  complex  conjugate  transpose 

tr,p(  ),<rmax(  )  trace,  spectral  radius,  largest  singular  value 

IlffWIloo  sup  ^mix  ] 

wglR 

n,  m,  mo ,  l  positive  integers 

A,B,C,K  nxn,  nxm,  lx  n,  mxl  matrices 

Bo, Co,  Do,  F  n  x  mo,  mo  x  n,  mo  x  mo,  mo  x  mo  matrices 

2.  Preliminaries 

In  this  section  we  establish  key  definitions  and  notational  conventions  that  simplify  the  expo¬ 
sition  in  later  sections.  We  begin  with  the  definitions  of  positive  real  and  bounded  real  transfer 
functions  [11,17]. 

In  this  paper  a  real-rational  matrix  function  is  a  matrix  whose  elements  are  rational  functions 
with  real  coefficients.  Furthermore,  a  transfer  function  ’is  a  real-rational  matrix  function  each  of 
whose  elements  is  proper,  i.e.,  finite  at  s  =  oo.  Finally,  a  stable  transfer  function  is  a  transfer 
function  each  of  whose  poles  is  in  the  open  left  half  plane.  Note  that  the  space  of  stable  transfer 
functions  is  denoted  in  [44]  by  RHoo,  i.e.,  the  real-rational  subset  of  Hqo. 

A  square  transfer  function  G(s)  is  called  positive  real  [17,  p.  216]  if  1)  all  elements  of  G(s)  are 
analytic  for  Refs]  >  0  and  2)  G(s)  +  G*(s)  is  nonnegative-definite  for  Re[s]  >  0.  A  square  transfer 
function  G(s)  is  called  strictly  positive  real  [2,14]  if  1)  all  elements  of  G(s)  are  analytic  for  Refs]  >  0 
and  2)  G(ju)  +  G*(ju)  is  positive  definite  for  real  u>.  Finally,  a  square  transfer  function  G(s)  is 
strongly  positive  real  if  it  is  strictly  positive  real  and  D  +  DT  >  0,  where  D  =  G( oo).  Note  that 
strongly  positive  real  implies  strictly  positive  real,  which  further  implies  positive  real.  Furthermore, 
we  note  that  if  a  transfer  function  is  strictly  positive  real,  then  the  system  is  stable  and  dissipative. 
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Next,  we  give  the  definition  of  bounded  real.  A  transfer  function  H(a)  is  bounded  real  [17]  if  and 
only  if  1)  all  elements  of  H(s)  are  analytic  for  Re[»]  >  0  and  2)  I  —  H(ju)H*(ju>)  is  nonnegative 
definite  for  real  u.  Equivalently,  2)  can  be  replaced  by  [17,  p.  307]  2')  I—  H(a)H‘(a)  is  nonnegative 
definite  for  Re[s]  >  0.  Alternatively,  a  transfer  function  H(a)  is  bounded  real  if  and  only  if  H(s)  is 
stable  and  satisfies  ||  (s) || oo  <  1. 


Next  we  establish  some  notation  involving  state  space  realizations  of  transfer  functions.  Let 


[44] 


G(«)  ~ 


A  I  B 


C  i  D 


(2.1) 


denote  a  state  space  realization  of  G(a),  that  is,  G(a)  =  C(al  -  A)~1B  +  D.  If  G(s)  is  square  and 
det  D  it  0,  then 


G-1(s)  ~ 


A  —  BD~1C  ,  BD 


-n 


-D~1C  D~l 

Finally,  if  Gx(s)  =  Ci(aT  -  A\)~lBi  +  Di  and  Gj(s)  =  Ca(8l  -  A^'^Bt  +  D 2,  then 


(2.2) 


’  Aa 

0 

b2 

Gx(s)Gj(s)  ~ 

BiCa 

Ax 

;  BaDa 

•  Di  C2 

Cx 

DiDa 

(2.3) 


3.  Riccati  Equation  Characterizations  of  Positive  Real 
and  Bounded  Real  Transfer  Functions 

In  this  section  we  provide  explicit  connections  between  positive  real  and  bounded  real  transfer 
functions  and  their  associated  state-space  realizations.  Furthermore,  we  give  Riccati  equation 
characterizations  of  their  resulting  state-space  realizations.  Finally,  we  draw  connections  with  the 
well-known  positive  real  lemma  [11,17,23]. 

We  begin  with  a  result  [17]  that  relates  bounded  real  transfer  functions  to  positive  real  transfer 
functions  via  the  Cayley  (bilinear)  transform.  Throughout  the  paper  7  denotes  a  positive  number. 

Lemma  3.1.  If  j~1JI(a)  is  an  mxm  bounded  real  transfer  function  with  det[Jm — jr_l/f(s)]  ^  0 
for  Re[a]  >  0,  then 

G(s)  ^  [/m  -  7-1ff(*)]-1[L»  +  'T1ff(«)]  (3.1) 

is  positive  real.  Conversely,  if  G(s)  is  an  m  x  m  positive  real  transfer  function  such  that  G(s)  is 
analytic  for  Re[s]  >  0,  then 

7 ''His)  ±  [G(s)  -  /m][G(s)  -  Imf 1  (3-2) 
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is  bounded  real. 


Proof.  Suppose  7-1.ff(s)  is  bounded  real.  Since  det[/m  -  7-1.ff(s)]  #  0  for  Re[s]  >  0,  it 
follows  that  G(s)  is  analytic  for  Re[s]  >  0.  Then  with  G(s)  defined  by  (3.1)  it  follows  that  'j~1H(s) 
satisfies  (3.2).  Thus,  we  obtain  for  Re[s]  >  0 

T'HW *  (a)  =  [G(a)  -  /m][G(s)  +  7m][G*(s)  +  -  Jm]  <  Im,  (3.3) 


which  implies 

[G(.)  +  /„]-‘[G-(.)  +  /„]-*  <  [G(<)  -  /J->[G-(.)  -/„]-■  (3.4) 

or,  equivalently, 

[<?*(«)  +  /m](G(s)  +  Im]  >  (G*(s)  -  7m][G(s)  -  7m]  (3.5) 

which  further  implies  that  G(s)+G*(s)  >  0  for  Re[s]  >  0.  Conversely,  suppose  G(s)  is  positive  real. 
Then,  since  G(s)  is  assumed  to  be  analytic  for  Re[s]  >  0,  it  is  easy  to  show  that  det  [G(s)  +  /m]  ^  0 
for  Re[s]  >  0.  Therefore,  y~1ff(a)  defined  by  (3.2)  is  analytic  for  Re[s]  >  0.  Then  with  7 ~1H(s) 
defined  by  (3.2)  it  follows  that  G(a)  satisfies  (3.1).  Next,  for  Re[s]  >  0  we  obtain 

G(*)+<r(s)  =  (/„l-7-1^(«)r1[^+'r-1Jff(«)]+(/m+7*1^W](/m-7-I^(®)]-1  > 0.  (3.6) 

Forming  [Im  -  7-1J7(s)](3.6)[/m  -  7~1R'*(s)]  yields 

[Im  +'j-1H(a)][Im  -  7 +  [Im  ~  7 +  [/m  +  7"1^*W]  >  0, 


which  implies  Im  -  7  2H(a)Hm[a)  >  0  for  Re[s]  >0.  □ 

Next,  we  use  the  results  of  Lemma  3.1  to  establish  connections  between  the  state  space  real¬ 
izations  of  positive  real  and  bounded  real  transfer  functions. 


Proposition  3.1.  If  G(a)  is  a  positive  real  transfer  function  with  minimal  realization 


G(s) 


A  B 
C  D  ’ 


(3.7) 


then  the  bounded  real  transfer  function  7  1H(a)  defined  by  (3.2)  has  a  minimal  realization 


7 -'His) 


'A  B 

C  D 


(3.8) 


4 


where 


A  =  A  —  B(Im  +  D)~1C, 

(3.9) 

£  4  vSs(7m  +  D)"1, 

(3.10) 

C  =  %ft(Im  +  D)~xC, 

(3.11) 

D  —  (D  —  7m)(Z?  +  7m)-1. 

(3.12) 

Conversely,  if  7  x77(s)  is  an  m  x  m  bounded  real  transfer  function  such  that  det[7m (s)]  ^  0 
for  Re[sJ  >  0  and  with  minimal  realization 


A-JL 

C  i  D 


> 


then  the  positive  real  transfer  function  G(s)  defined  (3.1)  has  a  minimal  realization 

A  B 

C  J5j» 

where 


(3.13) 


A  =  A  +  B(Im  —  D)~1C, 

(3.14) 

B  =  y/2B(Im  -  D)”1, 

(3.15) 

C  4  y/2(Im  -  D)~lC , 

(3.16) 

(3.17) 

Proof.  Given  (3.7)  it  follows  that  the  realizations  of  G(s)  -  Im  and  G(s)  +  Im  are  given  by 


[G(*)  ~  U] 


A  B 
C  D-Im  ’ 


[G(s)  +  7m] 


A__  B 
C  D  +  Im 


Now,  since  G(s)  is  positive  real,  it  follows  that  D  +  DT  >  0  which  further  implies  that  Im  +  D  is 
invertible.  Next,  using  (2.2)  we  have 


[(?(.)  + Jm]-» 


A  -  B(Im  +  D)~1C  Bilm  +  D)-1’ 
.  -{I„  +  D)-*C  {D  +  Im)~l  m 


Using  (2.3),  it  now  follows  that  7  1H(s)  =  [G(s)  -  7m][G(s)  +  7m]  1  has  a  nonminimal  realization 

r  A-B(Im  +  D)~'C  0  Bilrn  +  D)-1  1 

7 -B{Im  +  D)-*C  A  B(Im  +  D)~l 

■  (Im  ~  D)(Im  -f  D)~1C  C  (D  -  Im)(D  ~  Im)~l . 
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Next  it  follows  from  state-space  manipulations  that  7  1 
given  by 

\  A  —  B(Im  +  D)~1C 


H (s)  has  a  minimal  state-space  realization 
y/tBiln  +  D)-1  I 


7 -'His) 


L  n/2 {Im  +  D)-'C  ( D-Im)(D  +  Imy 

Furthermore,  Lemma  3.1  implies  that  /y~1ff(s)  is  bounded  real.  Finally,  the  converse  is  shown  in 
a  similar  fashion.  □ 


Having  established  connections  between  state-space  realizations  of  positive  real  and  bounded 
real  transfer  functions  we  proceed  in  the  spirit  of  recent  Hqo  results  [32-37]  to  establish  Riccati 
equation  characterizations  of  positive  real  systems. 


Theorem  3.1.  Let-ff(s) 
definite  matrix  <5  satisfying 


[A 

B 

C 

D 

where  trmtx(D)  <  7.  If  there  exists  annxn  nonnegative- 


0  =  AQ  +  QAr  +  7“2(B£>t  +  QCT)(Im  -  7 -iDDr)~1(BDT  +  QCr)r  +  BBr,  (3.18) 


then  (A,  B )  is  stabilizable  if  and  only  if 

A  is  asymptotically  stable.  (3-19) 


Furthermore,  in  this  case, 

II* tolloo  ^  7.  (3.20) 

Conversely,  if  A  is  asymptotically  stable  and  ||fl’(s)|j00  <  7,  then  there  exists  a  unique  nonnegative- 
definite  matrix  Q  satisfying  (3.18)  and  such,  that  the  eigenvalues  of  A+'j~3BDT(Im-^~2DDT)~1C 
+  7 -2QCT(Jm  —  7 _2DDt)-1C  lie  in  the  open  left  half  plane.  Furthermore,  Q  is  the  minimal 
solution  to  (3.18). 

Proof.  The  asymptotic  stability  of  A  follows  directly  from  Lyapunov  theory  while  (3.20) 
follows  from  algebraic  manipulation  of  (3.18);  for  details  see  [36].  The  converse  follows  from  the 
bounded  real  lemma  [17,  p.  308]  or  from  spectral  factor  theory  [31].  Finally,  the  proof  of  minimality 
is  given  in  [45].  □ 

Next,  we  utilize  a  transformation  that  converts  a  nonstrictly  proper  transfer  function  into  a 
strictly  proper  transfer  function  both  of  which  satisfy  the  same  Hoo  bound.  For  convenience  in 
stating  this  result  define  the  notation 

M  =  Im-1-7DDT,  N  =  Im  -  7"2£>t£. 
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Note  that  M  is  positive  definite  if  and  only  if  N  is  positive  definite. 


Proposition  3.2.  Let  H(s) 
is  asymptotically  stable  and 


A  B 
C  '  D 


and  assume  M  and  .V  are  positive  definite.  Then  A 


IlffMIloo  <  7 


(3.21) 


if  and  only  if  A‘  is  asymptotically  stable  and 


Ilff'WIU  <  i, 


where  if'(s) 


Aj  ,  Br 
C'  1  0 


and 


A1  =  A  +  't~2BDTM-1C, 
B'  =  BN ~1/2, 

C'  =  M~1/2C. 


(3.22) 


(3.23) 

(3.24) 

(3.25) 


Furthermore,  (3.18)  is  equivalent  to 

0  =  A'Q  +  QA/T  +  <y_aQC/TC,Q  +  B'B /T.  (3.26) 

Proof.  The  results  follow  from  Theorem  3.1  and  algebraic  manipulation.  For  details  see  [36]. 

□ 

Next,  using  Theorem  3.1  we  give  a  Riccati  equation  characterization  of  positive  real  trans¬ 
fer  functions.  To  do  this  we  use  (3.18)  to  imply  that  the  transfer  function  corresponding  to 

.  *  a  *  A  , 

[A,~fB,C ,~iD)  has  Hqo  norm  less  than  7.  By  Lemma  3.1  and  Proposition  3.1  the  resulting  Riccati 
equation,  i.e.,  (3.18)  with  ( A,B,C,D )  replaced  by  (A,~iB,C,iD),  implies  that  G[s)  ~ 
is  positive  real.  To  utilize  Theorem  3-1  we  require  that  <rta*x(,'lD)  <  7  or,  equivalently, 

Im  -  DD*  >  0.  (3.27) 

Now,  using  (.'.12),  this  is  equivalent  to 

Im  —  {D  -  Im)(D  +  Im)~1{D  +  Im)-r(D  -Im)T>  0.  (3.28) 

Since  {D  +  Im)  and  ( D  -  Im )  commute,  (3.28)  implies 

(Im  -  D)(Im  -  Dt)  <  (Im  t  D)(I  *  Dr),  (3.29) 


A  B 
C  D 
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which  implies  that 


D  +  Dt  >  0. 


(3.30) 


Thus,  we  restrict  our  attention  to  strongly  positive  real  systems. 

,  assume  detfD  +  7m]  ^  0,  define  A,B,C,D  by  (3.9)- 
(3.12),  and  assume  crm*x(D)  <  1.  If  there  exists  an  n  x  n  nonnegative-definite  matrix  Q  satisfying 

0  =  AQ  +  <2 At  +  ( BDt  +  <2 CT)(/m  -  DbT)-x(BbT  +  fiCT)T  +  BBr,  (3.31) 

A  A 

(A,  B)  is  stabilizable,  and 

det[/m  -  C(sln  -  A)~lB  -  D]  #  0  for  Refs]  >  0,  (3.32) 

then 

G(s)  is  strongly  positive  real.  (3.33) 

Conversely,  if  A  is  asymptotically  stable  and  G(a)  is  positive  real,  then  there  exists  a  unique 
nonnegative-definite  matrix  Q  satisfying  (3.31). 

Proof.  The  result  is  a  direct  consequence  of  Theorem  3.1,  Proposition  3.1  and  Lemma  3.1.  □ 

Remark  3.1.  Using  Proposition  3.2  we  can  represent  (3.31)  in  the  equivalent  form 

0  =  A'Q  +  QA'+  QC,tC'Q  +  B'B'r,  (3.34) 

where 

A'  =  A  -  B[Im  +  D)~XC  +  QB(Im  +  D)-1  DT (Im  -  DD T)~l(Im  +  D)~lC ,  (3.35) 

B'  =  y/2. B(Im  +  D -  brb)~1'2,  (3.36) 

6'  ±  V2(Im  -  bbT)-1/2(Im  +  D)~XC.  (3.37) 

Remark  3.2.  An  interesting  special  case  of  Theorem  3.2  is  the  case  D  =  Im.  Since  D  —  0 
(see  (3.12)),  (3.31)  or,  equivalently,  (3.34)  becomes 

0  =  {A-  i BC)Q  +  Q{A  -  \BC)T  +  \qC*CQ  +  ^BBT.  (3.39) 

it  it  4*  m 

Finally,  we  draw  connections  between  Theorem  3.2  and  the  well-known  positive  real  lemma 
used  to  characterize  positive  realness  in  the  state-space  setting  17 


Theorem  3.2.  Let  G(s)  ~ 
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A  B 


Lemma  3.2.  Let  G(s)  ~  I-—  — — j  be  an  m  x  m  transfer  function  with  minimal  realization 
(A,  B,C,  D).  Then  G(s)  is  positive  real  if  and  only  if  there  exist  matrices  Q  G  IRnXn,  L  €  IRnxp, 
and  W  G  IRmxp  with  Q  positive-definite  and  such  that 


0  =  AQ  +  QAt  +  LLr, 

(3.40) 

0  =  QCr  -B  +  LWt, 

(3.41) 

0  =  D  +  Dt  -  WWr. 

(3.42) 

This  form  of  the  positive-real  lemma  is  the  dual  of  that  given  in  [11,23],  and  the  derivation  is 
similarly  dual.  See  [12]  for  further  details  on  the  dual  positive  read  lemma. 

■  The  key  question  of  interest  here  is  the  relationship  between  Q  satisfying  (3.40)-(3.42)  and 
Q  given  by  (3.31).  To  answer  this  question,  we  invoke  the  assumption  that  D  +  DT  >  0  which, 
as  noted  earlier,  is  needed  for  the  existence  of  Q.  Thus,  once  again,  we  restrict  our  attention  to 
strongly  positive  real  transfer  functions.  In  this  case,  it  follows  from  (3.42)  that 


WWt  =  D+Dt.  (3.43) 

Now,  since  D  +  DT  >  0,  W  is  nonsingular  and  thus  (3.41)  implies 

L={B-QCt)W~t.  ,  (3.44) 

Using  (3.44)  it  follows  from  (3.40)  that 

0  =  AQ  +  QAt  +  {B-  QCt)W~tW~1(Bt  -  CQ)  (3.45) 

or,  since  (WW7)-'  =  W~TW~\ 

0  =  AQ  +  QAr  +  {B-  QCt){D  +  DT)~l{B  -  QCf.  (3.46) 


Thus,  we  have  shown  that  under  the  assumption  that  D  +  DT  >  0,  conditions  (3.40)-(3.42)  are 
equivalent  to  one  Riccati  equation  given  by  (3.46).  A  similar  result  for  the  dual  case  appears  in 

[17]- 

The  next  lemma  connects  the  two  Riccati  equations  (3.31)  and  (3.46). 

Proposition  3.3.  Assume  D  +  DT  >  0.  Then  the  Riccati  equation  (3.46)  is  identical  to  the 
Riccati  equation  (3.31),  or,  equivalently,  (3.34). 
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Proof.  Using  (3.46)  it  follows  that 


0  ={A-B(D  +  D^-'CjQ  +  Q[A  -  B(D  +  ZJT)-1C]T 

+  QCt{D  +  DT)~lCQ  +  B(D  +  Dt)~1Bt. 
The  result  now  follows  from  algebraic  manipulation  by  noting  that 


(3.47) 


(■ D  +  i?T)'1  =  2 (7m  +  D)-T[7m  -  (D  -  Im)(D  +  Im)-l{D  +  Im)~T(D  -  Jm)T]-1(/m  +  D)~\  □ 


Remark  3.3.  Note  that  in  the  case  D  =  Jm,  Proposition  3.3  can  readily  be  seen  by  comparing 
(3.39)  and  (3.46). 

4.  Robust  Stability  Problem  with  Positive  Real  Uncertainty 

In  this  section  we  state  the  robust  stability  problem  with  positive  real  uncertainty.  Consider 
the  uncertain  system 

x(t)  =  [A-B0F(Im  +  DoF^CoWt), 

(£) 

F+  Fr  >  0, 

when  the  inverse  of  7m  +  DoF  exists.  It  is  useful  to  note  that  (£)  can  be  viewed  as  a  strongly 
positive  real  system  ( A ,  Bo,C0,  D0)  in  a  negative  feedback  configuration  with  the  gain  F  (see  Figure 
4.1).  That  is, 


i(t)  =  Ax{t)  +  B0u(t),  (4.1) 

y(0  =  c0x(t)  +  D0u(t),  (4.2) 

with  negative  feedback 

u{t)  =  —Fy(t),  (4.3) 

where 

F  +  FT>0.  (4.4) 

Thus,  the  question  of  interest  is  the  stability  of  the  uncertain  system  (B)  with  positive  r eal  uncer¬ 
tainty  (4.4).  However,  before  we  proceed  with  this  question  we  give  a  lemma  on  the  existence  of 
(7m  +  DqF)~1  when  the  system  (A,  So,  Co,  Dq)  is  strongly  positive  real. 

Lemma  4.1.  Let  Dq,F  €  IRmoXm°)  and  assume  that  Dq  +  Dq  is  positive  definite  and  F  +  Ft 
is  nonnegative  definite.  Then 

det[7m  -I-  DoT’l  #  0.  (4.5) 
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Proof.  Since  Do  +  Dj  is  positive  definite,  it  follows  from  Lyapunov  stability  theory  that  -Do 
is  asymptotically  stable.  Hence  Do  is  nonsingular.  Furthermore,  it  follows  that  Dq1  4-  D^r  = 
Dq1  {Do  +  Dj)D0-T  is  positive  definite.  Since  F  +  FT  is  nonnegative  definite,  it  follows  that 
Dq1  +  F  +  (Do1  +  F)T  is  also  positive  definite.  Hence  -( Dq 1  +  F)  is  asymptotically  stable  and, 
consequently,  Dq1  +  F  is  nonsingular.  Thus 

det[/m  +  DoF]  =  (det  DoJde^Djj'1  +  F )  ^  0.  □ 

Next,  we  present  the  main  result  of  this  section  which  shows  that  the  uncertain  system  (D) 
is  robustly  stable  for  all  positive  real  uncertainty  of  the  form  (4.4).  For  the  statement  of  the  next 
result  we  define  the  uncertainty  set 

U  =  {AAeTRn*n:  AA  =  -B0F(Im  +  DoF)-1C0,  F  +  FT  >  0},  (4.6) 

where  Bo  G  IRnxmo,  Co  €  IRm°Xn,  and  Do  G  IRm°xm°  are  fixed  matrices  denoting  the  structure 
of  the  uncertainty  and  F  G  ntm°Xmo  jg  an  uncertain  matrix  (see  Figure  4.2). 

A  R  1 

Theorem  4.1.  Let  G(s)  ~  _  j  n°  ,  where  A  G  IRnXn  is  asymptotically  stable.  If  G(s) 

Oo  i  L/q 

is  strongly  positive  real,  then  A  +  A  A  is  asymptotically  stable  for  all  A  A  G  U.  Conversely,  if 
Do  +  Dq  >  0  and  A  +  A  A  is  asymptotically  stable  for  all  A  A  G  U,  then  G(s)  is  strongly  positive 
real. 

Proof.  As  shown  in  [2,14,18,21]  a  negative  feedback  configuration  consisting  of  a  positive  real 
transfer  function  and  a  strictly  positive  real  transfer  function  is  stable.  Under  the  assumption  that 
C(s)  is  strongly  positive  real,  the  dynamics  matrix  of  the  closed-loop  system,  which  has  the  form 
A  -  BoF(Im  +  D0F)-lCo,  is  asymptotically  stable.  Hence  A  -f-  AA  is  asymptotically  stable  for  all 
AA  G  U.  Conversely,  if  D0  +  Dq  >  0  and  A  +  AA  is  asymptotically  stable-for  all  AA  G  U,  then  it 
follows  from  the  definition  of  hyperstability  [21]  that  G(s)  is  strongly  positive  real.  □ 

Remark  4.1.  The  key  feature  of  the  uncertainty  set  U  is  that  the  uncertain  perturbation  AA 
involves  a  phase  constraint.  To  see  this  note  that  if  Do  +  Dq  >  0  and  F  +  FT  >  0,  then 

F(Im  +  DoF)"1  +  [F(/m  +  DoF)"1]1,  =  (/+  D0F)-T[F+  FT  +  FT(D0  +  Dj)F)(I+  DoF)-1  >  0. 

However,  the  term  F(/m  +  DoF)-1  is  bounded  in  magnitude  even  though  F  is  not.  For  example,  if 
F  is  a  scalar,  then  |F(l+D0F)-1|  <  1/Do-  Thus  the  uncertainty  set  U  incorporates  both  magnitude 
and  phase  constraints. 
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Remark  4.2.  Theorem  4.1  implies  that  robust  stability  of  the  uncertain  system  (17)  is  equiva¬ 
lent  to  a  positive  real  condition.  This  fact  can  be  compared  to  the  results  of  [37]  where  it  is  shown 
that  the  existence  of  a  fixed  Lyapunov  function  for  the  uncertain  system 


i(t)  =  (A  -t-  BoFCo)x(t), 

&m*x(F)  <  1, 


(£') 


is  equivalent  to  a  small  gain  condition.  However,  since  the  present  result  involves  a  phase  constraint 
not  present  a  small  gain  condition  (see  Remark  4.1)  one  should  expect  to  find  a  parameter  dependent 
Lyapunov  function  for  the  uncertain  systems  (i7)  rather  than  a  single  Lyapunov  function  as  in  [37]. 

5.  Robust  Controller  Synthesis  for  Positive  Real  Uncertainty 


In  this  section  we  state  the  Robust  Stabiliz ability  Problem  with  Positive  Real  Uncertainty.  The 
problem  involves  the  set  U  given  by  (4.6)  of  uncertain  perturbations  AA  of  the  nominal  (A,  S,C) 
system.  The  goal  of  the  robust  stability  problem  is  to  determine  a  static  output  feedback  controller 
that  stabilizes  the  plant  for  all  variations  in  U.  See  Figure  5.1. 


Robust  Stabilizability  Problem  with  Positive  Real  Uncertainty.  Determine  K  G 
IRmXt  such  that  the  closed-loop  system  consisting  of  the  nth-order  controlled  plant 


x(t)  =  (A  +  AA)x(t)  +  Bu[t),  t  G  [0,  oo), 


(5.1) 


measurements 

y(t)  =Cx(t), 


and  output  feedback  controller 

u(t)  =  Ky(t), 

is  asymptotically  stable  for  all  AA  €  U. 

For  each  uncertain  variation  AA  G  U,  the  closed-loop  system  can  be  written  as 

x(t)  =  (A -(- SJfC  +  ^A)x(t),  te[0,oo). 


(5.2) 


(5.3) 


(5.4) 


The  following  result  gives  necessary  and  sufficient  conditions  for  constructing  a  feedback  gain  K 
that  solves  the  Robust  Stabilizability  Problem  with  Positive  Real  Uncertainty.  For  the  statement 
of  this  result  define 

v  =  QCt(CQCt)-1C,  =  v. 
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Ro  =  (Do  +  D*)-\ 

for  arbitrary  Q  €  !RnXn  such  that  det  CQCT  ±  0,  and  let  Ri  and  i?2  be  arbitrary  real  n  x  n  and 
m  x  m  positive-definite  matrices. 

Theorem  5.1.  There  exists  K  €  IRmxi  that  solves  the  Robust  Stabilizability  Problem  with 
Positive  Real  Uncertainty  if  and  only  if  there  exist  n  x  n  positive-definite  matrices  Q,  P  satisfying 

0  =  (A  —  BR;lBrPv  -  BqRoC0)Q  +  Q(A  -  BR^BTPu  -  B0RoC0)T 

+  QCfRoCoQ  +  B0RoBf,  (5.5) 

0  =  (A  —  BqRqCq  +  QCq  RoCo)rP  +  P(A  —  BqRqCq  +  QCq  Po^o)  +  Pi 

-  PBP,1BTP  +  vJPBPf1BTP^±.  (5.6) 

Furthermore,  one  such  gain  K  is  given  by 

K  =  -R^1BTPQCT(CQCr)-1.  (5.7) 

Proof.  (Sufficiency).  Suppose  there  exist  n  x  n  positive-definite  matrices  Q,  P  satisfying  (5.5) 
and  (5.6).  Then,  with  K  given  by  (5.7),  it  follows  that  (5.5)  is  equivalent  to 

0  =  (A  +  BKC  -  B0RoC0)Q  +  Q(A  +  BKC  -  B0RoC0)T  +  QCJRoCqQ  +  B0RoBj,  (5.8) 

which  further  implies 

0  =  (A  +  BKC)Q  +  Q(A+BKC)r  +  (B0  -  QC?)T(D0  +  Dj^^Bo  -  QC0T)T.  (5.9) 
Furthermore,  (5.6)  is  equivalent  to 

0  =  (A  +  BKC  -  BoRqCq  +  QCf  RoC0)r  P  +  P(A+  BKC  -  B0RoC0  +  QCf  R0C0)  +  RX  +  KTRtK. 

(5.10) 

Note  that  (5.10)  is  an  auxiliary  equation  and  is  only  needed  for  computing  the  gain  K.  Furthermore, 
note  that  (5.9)  is  equivalent  to  (3.46),  or,  equivalently  (3.31).  It  now  follows  from  Theorem  3.2  that 
(A+BKC,  Bo, Co,  Do)  is  strongly  positive  real  which,  by  Theorem  4.1,  implies  that  A+BKC+AA 
is  asymptotically  stable  for  all  A  A  €  U. 

(Necessity).  It  follows  from  Theorem  3.2  and  Proposition  3.3  that  (A  +  BKC,  B0,C0,  D0)  is 
strongly  positive  real  if  and  only  if  there  exists  a  nonnegative-definite  solution  Q  to 

0  =  (A  +  BKC)Q  +  Q(A  +  BKC )T  +  (B0  -  QCj)(D0  -  Djy^Bo  -  QCj)T.  (5.11) 
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Now  it  follows  from  compactness  arguments  that  the  functional  J(K)  =  tr  Q(Ri  +  CT KT RjKC) 
must  have  a  global  minimum  on  the  set 

$  =  {if  €  IRmXi:  A  +  BKC  is  asymptotically  stable} 

under  the  assumption  that  BqBq  is  positive  definite.  Note  that  S  is  not  empty  by  the  assumption 
that  a  robustly  stabilizing  K  exists.  In  this  case  the  necessary  conditions  for  optimality  of  J(if), 
which  are  equivalent  to  the  existence  of  Q,  P  satisfying  (5.5)— (5.6),  must  have  a  solution.  □ 

Next,  we  specialize  Theorem  5.1  to  the  full  state  feedback  case.  When  the  full  state  is  available, 
i.e.,  C  =  the  projection  u  =  In  so  that  v±  =  0.  In  this  case  (5.7)  becomes 

K  -  -R^BTP  (5.12) 

and  (5.5),  (5.6)  specialize  to 

0  =  ( A-BRZlBrP-BoRoCo)Q+Q(A-BR;1BTP-B0RoCo)'T+QC2R0C0Q+BoRoBj , 

(5.13) 

0  =  (A  -  BqRqCq  +  QCfRoC0)rP  +  P(A  -  B0RoC0  +  QCJRqCq)  +  Rt  -  PBR^B^P. 

(5.14) 

It  is  interesting  to  note  that  even  in  the  full  state  feedback  case  the  result  involves  two  coupled 
Riccati  equations. 

A  salient  feature  of  (3.39)  is  the  fact  that  the  shift  —\BC  to  the  matrix  A  can  be  nonpositive. 

That  is,  -\BC  can  represent  a  left  shift  in  contrast  to  the  usual  a-shift,  which  is  a  uniform  open- 

loop  right  shift  used  to  place  the  closed-loop  poles  to  the  left  of  -a,  where  a  >  0[22] .  The  use  of 
a  left  shift  to  the  plant  dynamics  matrix  has  been  used  to  model  frequency  uncertainty  in  lightly 
damped  flexible  structures  [35-40].  Specifically,  consider  modal  dynamics  of  the  form 

A  =  block-diag f  Wl  ],...,  ~T)r  Wr  ]V  (5.15) 

\j-wx  —tii  [-wr  -T)r\J 

where  r\±  >  0  denotes  the  decay  rate  and  a/,-  denotes  modal  frequency.  Also  consider  uncertainty  of 
the  form 

r 

4A  =  5>A<,  (5.16) 

t=i 
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where  €  [-5, ,6,],  »'  =  1  are  real,  uncertain  parameters  with  given  bounds  6^,  and  the 

matrices  Ai  are  defined  by 


Ai 


block-diag (0, . . . ,  0, 


,0,  —  ,  0) , 


(5.17) 


where  the  matrix 


0  1 

-1  0 


corresponds  to  the  »th  diagonal  block  of  A.  The  skew  symmetric 
structure  of  Ai  accounts  for  uncertainty  in  the  ith  modal  frequency  w».  In  [10,  38-43]  the  Maximum 
Entropy  design  approach  is  predicated  upon  a  modified  covariance  (Lyapunov)  equation  of  the  form 


0  =  (A  +  S)Q  +  Q(A  +  S)T  +  ]T  SfAiQA?  +  V,  (5.18) 

i=l 


where  the  shift  5  is  defined  by 

*  »=1 

Note  that  5  has  the  form 

S  -  block-diag(— . . . , I2) 

so  that  S  effectively  shifts  each  mode  to  the  left  by  introducing  a  (fictitious)  augmentation  to 
the  open-loop  damping.  To  relate  (5.18)  to  (3.39),  consider  the  case  of  a  single  uncertain  modal 
frequency  by  setting  r  =  1.  Furthermore,  let 


*-*-«■*•  Al=[-i  oj- 

so  that  (with  B,C  replaced  by  B0,C0  in  (3.39))  -\B0C0  =  =  \S\A\  =  S.  The  remaining 

terms  5? A^QAj  +V  in  (5.18)  can  be  shown  to  play  a  role  similar  to  the  terms  QCTCQ  +  \BBt 
in  (3.39).  See  [46]  for  further  details.  Finally,  the  uncertain  perturbations  AA  given  by  (4.6)  have 
the  form 

AA  =  -S\F{h  +  DoF)-1.  (5.19) 

In  the  limiting  case  D0  — *  0,  setting  F  =  -jiAi  (so  that  F  +  Fr  >  0),  (5.19)  becomes 


AA  =  <Ti  Ai. 


Hence  U  given  by  (5.16)  can  be  used  to  capture  frequency  uncertainty  of  the  form  (5.16). 
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Abstract 

Nonlinear  controllers  offer  significant  advantages  over  linear  controllers  in  a  variety  of  cir¬ 
cumstances.  Hence  there  has  been  significant  interest  in  extending  the  linear-quadratic  synthesis 
methodology  to  nonlinear-nonquadratic  problems.  The  purpose  of  this  paper  is  to  review  the  cur¬ 
rent  status  of  such  efforts  and  to  present,  in  a  simplified  manner,  some  of  the  basic  ideas  underlying 
these  results.  In  particular,  we  focus  on  the  classic  paper  by  Bass  and  Webber  (1966)  and  show 
its  relationship  to  some  results  of  Speyer  (1976). 
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1.  Introduction 


Linear-quadratic  (LQ)  control  theory  has  been  extensively  developed  over  the  past  thirty  years. 
In  its  most  fundamental  form,  linear-quadratic  control  is  based  upon  the  following  assumptions: 

t)  the  plant  dynamics  and  measurement  equations  are  linear  in  both  the  state  and  control 
variables, 

it)  the  performance  measure  to  be  minimized  is  quadratic  in  the  state  and  control  variables, 

tit)  the  plant  disturbances  and  measurement  noise  are  additive  Gaussian  white  noise  or  Lt 
signals. 

In  addition  to  these  explicit  assumptions  the  following  implicit  assumptions  are  crucial: 

iv)  the  plant  model  is  completely  accurate, 

v)  the  state  and  control  variables  are  constrained  in  a  mean-square  or  Lj  sense. 

Under  these  assumptions  it  is  well  known  that  the  optimal  feedback  controller  is  linear  [l]. 

In  many  practical  situations,  however,  one  or  more  of  these  assumptions  may  be  violated.  For 
example,  the  plant  and  measurement  equations  may  be  nonlinear,  the  performance  measure  may 
be  nonquadratic,  the  disturbances  may  be  nongaussian  or  nonadditive,  the  plant  model  may  be 
uncertain,  or  state  variables  and  control  effort  may  be  limited  by  nonquadratic  constraints.  In 
such  cases  there  is  no  reason  to  expect  that  the  optimal  controller  iB  linear.  Rather,  it  should  be 
expected  that  nonlinear  controllers  will  have  better  performance  than  the  best  linear  controllers. 
For  example,  if  the  plant  model  is  nonlinear  then  nonlinear  controllers  can  be  used  to  account  for 
the  global  behavior  of  the  plant  [2-4].  Similarly,  gain-scheduled  controllers  designed  for  multiple 
plant  linearizations  constitute  a  widely  used  class  of  nonlinear  controllers  [5]. 

In  the  case  of  optimal  Hoo  performance  or  robust  stabilizability  in  the  presence  of  unstructured 
plant  uncertainty,  it  has  been  shown  [6—8]  that  nonlinear  controllers  offer  no  advantage  over  linear 
controllers.  However,  if  the  plant  uncertainty  is  structured  and  if  a  quadratic  Lyapunov  function 
is  assumed,  then  discontinuous  nonlinear  controllers  have  been  shown  to  offer  advantages  over 
linear  controllers  [9-14].  Continuous  approximations  to  the  discontinuous  controllers  of  [9,10]  have 
been  developed  in  [15,16].  Discontinuous  controllers  are  also  the  focus  of  variable  structure  control 
which  also  addresses  the  problem  of  plant  uncertainty  [17-20].  It  is  also  shown  in  [21]  that  nonlinear 
controllers  can  provide  improved  performance  in  a  neighborhood  of  the  worst  case  Hoo  disturbance 
attenuation. 
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Adaptive  controllers  can  be  viewed  as  nonlinear  controllers  that  operate  in  the  presence  of 
significant  plant  uncertainty.  Such  controllers  have  been  shown  to  stabilize  uncertain  systems  that 
cannot  be  stabilized  by  means  of  linear  controllers  [22-29]. 

With  regard  to  state  and  control  constraints,  one  of  the  most  common  nonlinearities  arising 
in  applications  is  actuator  saturation  [30-32].  Linear-quadratic  techniques,  however,  can,  at  best, 
only  impose  bounds  on  the  Lj  norms  of  the  state  and  control  variables.  Enforcing  constraints  of 
the  form  ||x(t)||  <  a  or  ||u(t)||  <  0  pointwise  in  time  requires  nonlinear  controllers. 

In  view  of  the  advantages  of  nonlinear  controllers  over  linear  controllers,  it  is  not  surprising  that 
significant  effort  has  been  devoted  to  developing  a  theory  of  optimal  nonlinear  regulation  [33-70]. 
The  goal  of  the  present  paper  is  to  provide  a  simplified  framework  for  optimal  nonlinear  regulation 
in  terms  of  nonquadratic  cost  functionals.  In  accordance  with  practical  motivation,  we  restrict 
our  attention  to  time-invariant  systems  on  the  infinite  interval.  In  this  case  asymptotic  stability  is 
guaranteed  by  means  of  a  Lyapunov  function  for  the  closed-loop  system.  This  Lyapunov  function 
is  given  as  the  solution  to  a  steady-state  form  of  the  Hamilton- Jacobi- Bellman  equation. 

In  future  research,  we  intend  to  reverse  the  situation  somewhat  by  fixing  the  structure  of 
the  Lyapunov  function,  cost  functional,  and  feedback  law  prior  to  optimization.  In  this  case  the 
structure  of  the  Lyapunov  function  can  be  viewed  as  providing  the  framework  for  controller  synthesis 
by  guaranteeing  local  or  global  asymptotic  stability  for  a  class  of  feedback  controllers.  The  actual 
controller  chosen  for  implementation  can  thus  be  the  member  of  this  candidate  class  that  minimizes 
the  given  performance  functional.  In  LQG  theory,  for  example,  the  Lyapunov  function  is  the 
familiar  quadratic  functional  V{x)  =  xTPx,  while  the  gains  for  the  linear  feedback  control  are 
chosen  to  minimize  a  quadratic  performance  functional.  In  summary,  then,  Lyapunov  function 
theory  provides  the  framework,  while  optimization  fixes  the  gains. 

2.  Nonquadratic  Cost  Evaluation 

In  this  section  we  investigate  the  role  of  Lyapunov  functions  in  evaluating  nonquadratic  cost 
functionals.  To  expand  upon  the  linear-quadratic  case,  we  consider  the  problem  of  evaluating  a 
nonquadratic  cost  functional  depending  upon  a  nonlinear  differential  equation.  It  turns  out  that 
the  cost  functional  can  be  evaluated  in  closed  form  so  long  as  the  cost  functional  is  related  in  a 
specific  way  to  an  underlying  Lyapunov  function.  The  basis  for  the  following  development  is  the 
paper  [36]  by  Bass  and  Webber.  Note  that  the  results  of  this  section  make  no  explicit  reference  to 
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control. 


In  accordance  with  practical  motivations,  we  restrict  our  attention  to  time-invariant  systems  on 
the  infinite  horizon.  Furthermore,  for  simplicity  we  shall  define  all  functions  globally  and  assume 
that  existence  and  uniqueness  properties  of  the  given  differential  equations  are  satisfied. 

For  the  following  result,  let  / :  IR”  — *  IR”  and  L:  IRn  — +  IR.  We  assume  /( 0)  =  0.  Let  (') 
denote  derivative. 


Lemma  2.1.  Consider  the  system 

±(t)  =  f(x(t)),  z(0)  =  x0,  t>  0,  (2.1) 

with  performance  functional 

j(x o)  =  f°°  L(x(t))dt.  (2.2) 

J  o 

Furthermore,  assume  there  exists  a  C1  function  V :  IRn  — *  IR  such  that 

V  (0)  ss  0,  (2.3) 

V(x)>0,  x€lR",  x#0,  (2.4) 

V'(z)f(x)<  0,  x€  JR”,  x#0,  (2.5) 

L(x)  +  V'[x)f(x)  =  0,  *  €  IR”.  (2.6) 

Then  x(t)  =  0,  t  >  0,  is  a  globally  asymptotically  stable  solution  to  (2.1)  with  Xo  =  0.  Furthermore, 

J(x0)  =  V(*0),  x0  e  IR".  (2.7) 

Proof.  Let  x(t),  t  >  0,  satisfy  (2.1).  Then 

V(*(l))  =  Jt  V(*(.‘))  =  V'(»(i))/(l(l)).  ‘  >  0.  (2.8) 

Hence  it  follows  from  (2.5)  that 

V{x{t))  <0,  t  >  0,  x(t)  #  0.  (2.9) 

Thus,  by  (2.3),  (2.4),  and  (2.9)  it  follows  that  V(-)  is  a  Lyapunov  function  for  (2.1)  and  thus 
x(t)  -♦  0  as  t  —*  oo  for  all  initial  conditions  xo-  This  proves  global  asymptotic  stability  of  the 
solution  x(t)  =0,  t  >  0.  Now  (2.8)  implies  that 

0  =  -V(x(*))  +  V'{x{t))f(x(t)) 
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and  hence,  by  (2.6), 


£(*(<))  =  -v  WO)  +  £(*(0)  +  v'(x(0)/(i(()) 

=  -v(*(i)). 

Integrating  over  [0,  t)  yields 

/*  £(*(*)  )ds  =  -  V(x(t))  +  V  (*0). 

Jo 

Now  letting  t  —*  oo  and  noting  that  V(x(t))  — *  0,  yields  (2.7).  □ 

The  main  feature  of  Lemma  2.1  is  the  role  played  by  the  Lyapunov  function  V  (z)  both  in 
guaranteeing  stability  and  in  evaluating  the  functional  J(x 0).  Let  us  illustrate  this  result  with  a 
familiar  example.  Consider  the  linear  system 


with  cost  functional 


x  =  Ax,  x(0)  =  x0, 


J(xo)  =  /  xrRx  dt, 
Jo 


(2.10) 


(2.11) 


where  R  €  IRnxn  is  positive-definite.  If  A  is  asymptotically  stable  then  there  exists  a  positive- 
definite  matrix  P  €  IRnxn  satisfying 


0  =  ArP  +  PA+R. 


(2.12) 


Now  define 


V(x)  =  xTPx,  (2.13) 

which  satisfies  (2.3)  and  (2.4).  rthermore,  with  f(x)  =  Ax  and  L(x)  =  xrRx  it  follows  that 

V'(x)f(x)  =  2  xtPAx  =  xT  (AT  P  +  PA)x  =  -x^Rx  =  -L(x), 

which  verifies  (2.5)  and  (2.6).  Hence 

J(x0)  =  xfPx0, 

which  is  a  familiar  result  from  linear-quadratic  theory. 

Remark  2.1.  Note  that  if  (2.6)  is  valid,  then  (2.5)  is  equivalent  to 

L( x)  >0,  i  €  IRn,  x  #  0. 


(2.14) 
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More  generally,  assume  A  is  asymptotically  stable,  let  P  be  given  by  (2.12),  and  consider  the 

case 

L(x)  =  xrRx  +  h[x),  (2-15) 

f{x)  =  Ax  +  N(x),  (2.16) 

V  (i)  =  *  TPx  +  g(x),  (2.17) 

where  h(-)  and  <?(•)  are  nonquadratic  and  N(-)  is  nonlinear.  To  satisfy  (2.6)  we  require  that 

xrRx  +  h(x)  +  [2 xTP  +  5'(i)][Ai  4-  N(x)]  =  0.  (2-18) 

Our  goal  is  to  study  (2.18)  under  a  variety  of  choices  for  h(-),g(-),  and  N(-).  For  convenience, 
rewrite  (2.18)  as 

x  t(AtP  +  PA  +  R)  x  +  g'{x)Ax  +  h(x)  +  g'(x)N(x)  =  0.  (2-19) 

If  A  is  asymptotically  stable,  then  we  can  choose  P  to  satisfy  (2.12)  as  in  the  linear-quadratic  case. 
Next,  suppose  N(x)  =  0.  Then  (2.19)  is  satisfied  if 

g,(x)Ax  +  h(x)  =  0.  (2.20) 

The  following  lemma,  which  is  quoted  in  [36],  will  be  useful  for  satisfying  (2.20). 

Lemma  2.2.  Let  A  €  IRnXn  be  asymptotically  stable  and  let  hk:  IRn  — ►  IR  be  a  nonnegative- 
definite  homogeneous  A:-form  (A:  even).  Then  there  exists  a  unique  nonnegative-definite  homoge¬ 
neous  fc-form  gk’-  IR"  — *  IR  such 

gk(x)Ax  +  hk(x)  =  0,  *  €  IR".  (2-21) 

Proof.  The  result  can  be  shown  using  the  Kronecker  product  representation  of  multilinear 
functions.  Details  will  appear  in  an  expanded  version  of  this  paper.  □ 

Suppose  now  that  h(x)  is  of  the  form 

*(*)  =  Yi  (2-22) 

vs2 

where,  for  u  =  2,...,r,  IR”  — *  IR  is  a  nonnegative-definite  homogeneous  2i/-form.  Now, 

using  Lemma  2.2,  let  IRn  — ♦  IR  be  the  nonnegative-definite  homogeneous  2i/-form  satisfying 

g'2u{x)Ax  +  h2u(x)  =  0,  x  €  IRn,  *'  =  2,...,r,  (2.23) 
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and  define 


(2.24) 


Since 


$(*)  =  £  *»"(*)• 


$'(*)  =  52  **»(*), 


summing  (2.23)  over  u  yields  (2.20).  Since  (2.6)  is  now  satisfied,  (2.9)  implies  that 


J(x0)  =  xJPxq  +g(x0). 


(2.25) 


As  another  illustration  of  condition  (2.18),  suppose  that  V(z)  is  constrained  to  be  of  the  form 


V(x)  =  xTPx  +  -(xTAfx)a, 

m 


(2.26) 


where  P  satisfies  (2.12)  and  M  is  an  n  x  n  symmetric  matrix.  Then,  with  N(x)  =  0  (2.18)  yields 


h(x)  =  —(xtMx)xt(AtM+  MA)x. 


(2.27) 


If  R  is  an  n  x  n  symmetric  matrix  and  M  is  given  by 


0  =  ArM  +  MA  +  R, 


(2.28) 


then  h(x)  satisfying  (2.18)  is  of  the  form 


h(x)  =  (xr  Mx)(xr  Rx). 


(2.29) 


Thus,  if  V(x)  is  of  the  form  (2.26),  then,  by  utilizing  (2.28),  condition  (2.18)  is  satisfied  if  L(x)  has 


the  form 


L{x)  =  xtRx  +  (xTMi)(xTfix). 


(2.30) 


In  the  next  section  we  apply  Lemma  2.1  to  the  problem  of  optimal  nonlinear  feedback  control. 
The  relation  (2.18)  shall  play  a  key  role  with  greater  complexity  arising  from  the  fact  that  the 
nonlinear  dynamics  term  N(x)  will  be  nonzero. 

3.  Optimal  Control 

In  this  section,  we  extend  the  development  of  Section  2  to  obtain  a  characterization  of  optimal 
feedback  controllers.  These  conditions  are  essentially  a  specialization  of  the  Hamilton-Jacobi- 
Bellman  (HJB)  conditions  for  the  time-invariant,  infinite-horizon  case.  For  this  problem  the  HJB 
partial  differential  equation  reduces  to  a  purely  algebraic  relationship. 
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We  begin  with  a  notion  of  optimality  involving  only  a  Lyapunov  function.  Hence  let  / :  IR"  x 
IRm  —*  ntn,  assume  that  /(0,0)  =  0,  and  consider  the  controlled  system 

*W  =  /(*(*)»  “M)»  *(0)  =  *o,  t>  0.  (3.1) 

The  control  u(-)  in  (3.1)  is  restricted  to  the  class  of  admissible  controls  consisting  of  measurable 
functions  u:  (0,  oo)  — ►  IRm  such  that 

u(t)  e  n,  t  >  0,  (3.2) 

where  D  C  IRm  is  given. 

Definition  3.1.  Let  V :  IRn  — » IR  be  a  C1  function.  The  function  <f> :  IRn  -*fl  is  optimal 
with  respect  to  V  if 

V'(*)/(*^(*))<V'(x)/(*,u),  x€TRn,  u€f2.  (3.3) 

Note  that  if  V  satisfies  (2.3),  (2.4),  and 

V'(x)/(x,*(x))  <  0,  *€lRn,  x#0,  (3.4) 

then  V  is  a  Lyapunov  function  for  the  closed-loop  system 

x(t)  =  f(x(t),j(x(t)),  x(0)  =  x0,  t>  0,  (3.5) 

that  is,  system  (3.1)  with  u(t)  =  ^(z(t)),  t  >  0.  The  inequality  (3.3)  thus  characterizes  feedback 
controllers  that  optimize  the  decay  rate  of  the  closed-loop  system  as  measured  by  the  Lyapunov 
derivative. 

To  illustrate  Definition  3.1,  consider  the  linear  system 

i(t)  =  Ax(t )  +  Bu{t),  x(0)  =  xQt  t>  0,  (3.6) 

where  A  is  asymptotically  stable  and  the  control  constraint  set  Cl  is  characterized  by 

t>0,  t  =  l,...,m,  (3.7) 

for  given  positive  constants  ai, . . .  ,am.  Define  the  quadratic  Lyapunov  function 

V(x)  =  xrPx,  x  €  IRn,  (3.8) 
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where  P  is  the  unique  positive-definite  solution  to 

0  =  AtP  +  PA  +  R  (3.9) 

for  an  arbitrary  n  x  n  positive-definite  matrix  R.  Then 

V'(x)/(x,u)  = -xTRx  +  2(BTPx)Tu,  xeUT,  u  €  IRm.  (3.10) 

It  is  easy  to  see  that  a  feedback  control  <j>  —  (<j>i , . . . ,  ^m)T:  IRn  — »  f)  that  is  optimal  with  respect 
to  V  is  given  by 

4>i{x)  =  -a»sgn(BTPx),,  *  =  l,...,m.  (3.11) 

Note  that  if  (BT Px)i  =  0,  then  V'/(x)/(x,u)  is  independent  of  u^.  Hence  the  value  of  <f>i(x)  has  no 
effect  on  (3.3).  If  in  place  of  (3.7)  we  constrain  u(t)  by 

||«MII  <  a,  *  >  0,  (3.12) 

where  |j  •  ||  denotes  the  Euclidean  norm  and  a  >  0,  then  a  feedback  function  <f>  that  is  optimal  with 
respect  to  V  is  given  by 

*(l)  =  p*feji BTpx'  *  fiTpi  *  °'  (313) 

with  ^(x)  arbitrary  if  BTPx  =  0. 

We  now  turn  to  the  problem  of  characterizing  feedback  controllers  that  minimize  a  performance 
functional.  Let  L:  IRn  x  IRm  — » IR  and,  for  p  €  IRm,  define 

H(x,p,u)  =  L(x,u)+pT/(x,u). 

Theorem  3.1.  Consider  the  controlled  system  (3.1)  with  performance  functional 


J(xo,u(-))  =  L(x(t),u(t))dt.  (3.14) 

Jo 

Assume  that  there  exist  a  C1  function  V:  IRn  — ►  IR  and  a  function  4>‘  IRn  — *  fi  such  that 

^(0)  =  0,  (3.15) 

V (x)  >  0,  x  e  IRn,  x#0,  (3.16) 

^(0)  =  0,  (3.17) 

V"(x)/(x,*(x))<0,  xeIR",  x#0,  (3.18) 

tf(x,V'T(x),*(x))  =  0,  x  €  IRn,  (3.19) 

H(x,V'T(x),u)  >  0,  xenr,  u  €  I?-  (3.20) 
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Then,  with  the  feedback  control  «(•)  =  <f>[x(’)),  the  solution  x(t)  =  0,  t  >  0,  of  the  closed-loop 
system  is  asymptotically  stable,  and 

J(*o,*(x(.)))  =  V(x0).  '  (3.21) 

Furthermore,  the  feedback  control  «(•)  =  ^>(x(-))  minimizes  J(xo,u(-))  in  the  sense  that 

J(x0>(x(-)))=  min  J(x0, «(•)),  (3.22) 

u(-)€a(*0) 

where 

U(x o)  =  {«(•):  u(*)  is  admissible  and  x(*)  given  by  (3.1)  satisfies  lim  V(x(t))  =  0). 

t— *oo 

Proof.  Global  asymptotic  stability  and  (3.21)  are  obtained  by  using  (3.15)-(3.19)  and  applying 
Lemma  2.1  to  the  closed-loop  system  (3.5).  To  prove  (3.22),  let  tt(-)  &  U(x 0)  and  let  x(-)  be  the 
solution  to  (3.1).  Then  it  follows  that 

V(x(t))  =  V#(x(t))/(x(t),u(t)) 

or 

0=-V(*(t))  +  V'(x(t))/(*(t),«(t)). 

Hence 

£(*(<).“(<))  =  -V"(i(l))  -T  i(i(l),u(l))  +  v'(l(t))/(l(l),<*(0) 

=  -v(z(t))+  H(i(l),V'T(x(t)),«(i)). 

Now  using  the  fact  that  u(-)  €  U{x o)  along  with  (3.20)  and  (3.21),  we  obtain 

■'(*»,«(•))  =  fl-vW  i))  +  ^(*W.v,T(x(i)),«(t))ldi 

Jo 

=  -  lim  V(x(t))  +  V(x0)+  f  flr(*(0iv’'T(*W)»«W)' d* 

t-»oo  Jq 

=  V(*0)  +  r  F(*(t),V'T(x(*)),u(0)dl 
Jo 

>  V(x0) 

=  «7(*o,  <£(*(•))> 

which  yields  (3.22).  □ 

The  principal  feature  of  Theorem  3.1  is  that  the  optimal  control  law  u  =  <f>(x)  is  a  feedback 
controller.  Furthermore,  this  control  is  an  optimal  stabilizing  control  independent  of  the  initial 
condition  xq. 
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Remark  3.1.  If  (3.19)  and  (3.20)  are  satisfied,  then  it  follows  that 


L(x,<f>(x))  +  V'(x)f(x,<j>(x))  <  L(x,u)  +  V'(x)/(x,u),  x  €  IR",  u  €  fi.  (3.23) 

If  L(x,  u)  is  independent  of  u,  then  (3.23)  is  equivalent  to 

v\x)f  {*><!>{*))  {*,*)>  *€lRn,  uen,  (3.24) 

which  is  precisely  condition  (3.3).  Thus,  in  this  case  conditions  (3.19)  and  (3.20)  imply  that  the 
feedback  control  u(-)  =  ^>(x(-))  is  optimal  with  respect  to  V . 

Now  let  us  illustrate  Theorem  3.1  with  some  examples.  We  begin  with  the  simplest  case, 
namely,  the  linear-quadratic  regulator.  Hence  consider  the  controlled  system 

x  =  Ax  -f  -Btt,  x(0)  =  x0,  t  >  0,  (3.25) 

with  performance  functional 

J(x0,  «(•))=  f  [xTf?iX  +  uTR2«]dt,  (3.26) 

Jo 

where  x  €  IRn,  u  €  IRm,  and  where  Ri  and  l?j  are  positive  definite.  Thus,  L(x,u)  has  the  form 

L(x,  u )  =  xrRii  +  uTi?2tt-  (3.27) 

Furthermore,  assume  that  V (x)  is  quadratic,  that  is, 

V(x)  =  xTPx,  (3.28) 

where  P  is  positive  definite.  Consequently,  we  have 

H(x,ViT(x),u)  =  xr(ATP  +  PA  +  Ri)x  -(-  2xT PBu  +  ur R2U. 

If  (3.19)  and  (3.20)  are  satisfied,  then  u  =  <j>{x)  must  minimize  H(x,VfT{x),u).  Hence,  setting 

jZH{x,V't{x),u)  =  0  (3.29) 

yields 

u  =  ^(x)  =  -Rj-1HtPx.  (3.30) 

To  check  (3.19)  we  note  that 

H{x, V'(x),  ^(x))  =  xt(AtP  +  PA+Rl-  PSP)x , 
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where  5  =  BBj  1BT.  Thus,  (3.19)  holds  if  P  is  chosen  to  satisfy 


0  =  AtP  +  PA+Rx-  PSP. 

(3.31) 

The  closed-loop  system  has  the  form 

x  =  (A  -  SP)x,  x(0)  =  xo,  t  >  0. 

(3.32) 

Writing  (3.31)  as 

0  =  (A  -  SP)TP  +  P(A  -  SP)  +  Rlt 

(3.33) 

it  follows  that  A  — 

SP  is  asymptotically  stable.  Finally,  it  follows  using  (3.31)  that 

H(x,V,T(x),u)  =  (u  +  Bf 1  BTPx)  TB2  (u  +  Bf 1  BTPx) 

(3.34) 

so  that  (3.20)  is  satisfied.  In  summary,  the  solution  u  =  4>{x)  to  the  linear-quadratic  problem  is 
given  by  (3.30)  where  P  is  the  positive-definite  solution  to  (3.31). 

Next,  we  consider  the  case  of  a  nonquadratic  cost  and  nonlinear  feedback  control.  Hence  assume 
that  L(x,u),  /(x,u),  and  V(x)  are  of  the  form 


L(x,u)  =  xtBix  +  h(x)  +  uTBa  ti. 

(3.35) 

/(x,tx)  =  Ax  +  Bu, 

(3.36) 

V(x)  =  xTPx  +  g(x). 

(3.37) 

With  this  notation  we  have 

H(x,V'T(x),u)  =  xT  Rix  +  h(x)  +  utBju  +  [2xTP  4-  g'(x)][Ax  +  Bu]. 

Again  setting  ^H(x,V'T(x),u)  =  0  we  obtain 

u  =  4>(x)  =  -Rj1BrPx  -  if^BVfx).  (3.38) 

This  yields 

B(x, V'T (x) , ^(x))  =  xT{ATP+PA+Rl-PSP)x+h{x)+g'(x){A-SP)x-lg'{x)Sg,r{x).  (3.39) 

To  satisfy  (3.19),  let  P  satisfy  (3.31)  and  require  that 

h(x)  +  g\x){A-SP)x-  ?-g'{x)Sg'T{x)  =  0,  x  €  IR".  (3.40) 

4 
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With  (3.31)  and  (3.40)  it  follows  that  (3.19)  is  satisfied.  Furthermore,  it  can  be  shown  that 

if(z,F/T(x),u)  =  [u  -  ^(x)]Tf?j[u  -  <£(x)]  (3.41) 

so  that  (3.20)  holds. 

Returning  to  (3.40),  let  us  consider  the  approach  of  [36].  Suppose  that  for  v  =  2, . . . ,  r,  h2„(z)  is 
a  given  nonnegative-definite  homogeneous  2i/-form.  Since  A-  SP  is  asymptotically  stable,  Lemma 
2.2  implies  that,  for  v  =  2, . . .  ,r,  there  exists  a  nonnegative-definite  homogeneous  2i/-form  g2U(x) 
satisfying 

h2u{x)  +  g'iv{x){A- SP)x  =  0,  i/  =  2,...,r.  (3.42) 

Then  (3.40)  is  satisfied  with  h(x)  and  g(x)  defined  by 

M*)  =  h**ix)  +  ^5,(*)^,t(z)>  (3.43) 

v=2 

T 

*(*)  =  S***^*)-  (3-44) 

As  discussed  in  [36],  the  term  jg,(x)Sg,T(x)  appearing  in  h(x)  in  (3.43)  is  somewhat  artificial  in 
the  sense  that  it  cannot  be  specified  arbitrarily.  It  is  interesting  to  note,  however,  that  with  h(z) 
given  by  (3.43),  the  performance  functional  has  the  form 

roo  r 

J(z0,u(-))  =  /  [xT  Rix+y^hivfa)  +  uT  R}\i  + <j>JjL[x)R2<f>NL{x)]dt.  (3.45) 
J° 

In  (3.45),  4>nl{x)  is  the  nonlinear  part  of  the  optimal  feedback  control,  that  is, 

<t>{x)  —  4>l{x)  +  <f>NL(x),  (3.46) 

where 

<MI)  =  -RZ1BTPx,  4>nl{x)  =  ~R;1BTg'T(x).  (3.47) 

As  another  example,  suppose  we  require  that  V(x)  be  of  the  form 

V  (x)  =  xT  Px  +  |  (xT  M  x)  2 ,  (3 .48) 

which  corresponds  to  (3.37)  with  g(x)  =  |(xTMx)}.  Thus  (3.38)  specializes  to 


u  =  ^(x)  =  -Ry1  BT  Px  —  R2  1Bt(xtMx)Mx 


(3.49) 


and  (3.40)  becomes 

h(x)  -|-  (xTMx)xr[(A  -  SP)rM  +  M(A  -  SP)]x  -  {xTMx)ixr  MSMx  =  0.  (3.50) 

To  satisfy  (3.50),  let  R  be  an  arbitrary  n  x  n  symmetric  matrix  and,  since  A-  SP  is  asymptotically 
stable,  let  M  be  the  unique  symmetric  solution  to 

( A-SP)TM  +  M(A-SP)  +  R  =  0 .  (3.51) 

Now  (3.50)  is  satisfied  with 

h(x)  =  (xT  Mx)(xT  Rx)  +  (itMi),itMSMi  (3.52) 

A  stochastic  version  of  this  problem  was  treated  in  [49].  In  [49]  the  matrix  R  is  chosen  to  be 
Ri  +  MSM  so  that  (3.49)  becomes  a  Riccati  equation 

(A  -  SP)tM  +  M(A  -  SP)  +  Ri  +  MSM  =  0.  (3.53) 

See  equation  (23)  of  [49]  setting  W2  =  0. 
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